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ABSTRACT 


Aircraft parameter estimation is the process of extracting 
numerical values for aerodynamic stability and control derivatives 
from flight-test time history data. This process can be used as a 
verification or validation tool for results obtained from wind- 
tunnel testing or through computational analysis, and can obtain or 
improve estimations of dynamic derivatives. 

This study implements the MATLAB Personal Computer (PC) based 
maximum likelihood estimation routine for aircraft longitudinal and 
lateral-directional derivatives. The parameter estimation was 
first accomplished on generated simulated data, with and without 
noise. The noise consisted of measurement and state noise which 
used the Dryden Gust Model. Secondly, two actual longitudinal 
flight-test maneuvers are analyzed for the F-14A and the T-37 
aircraft. Additionally, the simulated portion of this study can be 
an excellent instructional aid in Flight Dynamics and Flight Test 
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I. INTRODUCTION 


Aircraft parameter estimation is the process of extracting 
numerical values for aerodynamic stability and control 
derivatives from flight-test time history data. Aircraft 
flight tests designed for this purpose are generally motivated 
by one or a combination of the following objectives: 

1. The desire to correlate flight test parameter estimates 
with wind tunnel data and analytic results. 


2. The desire to more accurately refine parameter estimates 
for purposes of control system analysis and design. 


3. The desire to achieve an accurate aircraft math model for 
use in high fidelity flight simulators. 

An early and continued use of parameter estimation, as stated 

above, is in the validation of wind tunnel and analytic 

results. However, due to the continuing advances in aircraft 

design and performance capabilities, the ability to accurately 

extrapolate wind-tunnel test results is diminishing and a 


greater emphasis is being placed on flight test results. 


feef.l:p.2] 
Comprehensive wind-tunnel testing, combined with 
analytic analysis, can give reasonable estimates of an 


aircraft’s aerodynamic derivatives, but there are potential 
sources for inaccurate predictions: the matching of "scaled" 


wind-tunnel tests with expected flight conditions is 


difficult, with Reynolds number differences often being the 
standard explanation for discrepancies. Reliable and accurate 
dynamic wind tunnel test results are extremely difficult and 
expensive to accomplish. Support systems (stings, etc.) have 
become an issue as to the extent the data, especially drag, 
are affected. [Ref. 2:p.3] Additionally, the ever present 
time and money constraints often necessitate shortcuts in not 
only wind-tunnel testing but in flight testing as well. It 
seems wise, therefore, to use flight-test data, at the very 
least, as a verification tool of aircraft stability and 
control derivatives for even the most simple configurations. 

Currently at the Naval Postgraduate School (NPS), in the 
Department of Aeronautics and Astronautics, research is being 
done uSing remotely piloted aircraft as research testbeds. 
One testbed in use has been a half-scale Pioneer Unmanned Air 
Vehicle (UAV). 

The full-scale Pioneer is operational in the U.S. Navy and 
Marine Corps and saw extensive action in the recent war with 
Iraq. The small size of the Pioneer or essentially any 
tactical UAV allows it to operate close to and in some cases 
behind enemy lines, extending the "eyes" of battlefield 
commanders. Its missions include gun fire spotting, real time 
enemy surveillance, bomb damage assessment, target designation 
and an array of intelligence collection techniques. 

The relatively low cost, small size, reduced risk and 


inherent flexibility of UAV’s, such as the half-scale Pioneer, 


have allowed the department to become actively involved in 
assessing their flight characteristics. The Pacific Missile 
Test Center (PMTC) at Pt. Mugu California is a development and 
testing facility for the U.S. Navy. Current activity at PMTC 
includes developmental work in conjunction with the Pioneer 
UAV. In development by the Target Simulation Lab at PMTC is 
a flight training simulator to be used by Pioneer operators 
for initial training and proficiency flights. Results from 
thesis work of two former NPS students, USMC Capts. Daniel 
Lyons and Robert Bray, have been supplied to the Lab [Refs. 3 
and 4]. These results comprised various aerodynamic 
parameters obtained from two different approaches, a numerical 
method (low-order panel technique) and wind tunnel tests of a 
0.4-scale model. The aerodynamic data supplied to PMTC are 
being used in their math model for the simulator. 

The goal of the present study is to incorporate a personal 
computer (PC) parameter estimation capability into the ongoing 
NPS flight research. This application will give the flight 
test program an added dimension: the ability to compare data 
from wind tunnel and analytic analyses with flight test 
results. Additionally, the adapted program can be 
incorporated into flight test and dynamic stability and 
control courses as a valuable teaching aid. 

In the near term, interest in the Pioneer parameter 
estimation results has been expressed by PMTC in hopes of 


achieving a more realistic training simulator for the 


operators. It is hoped that full scale time history data can 
be obtained from PMTC. Future work in this area includes 
completion of the Pioneer flight research and comparison of 
the derivative results obtained from time history parameter 
estimation with those obtained in References 3 and 4. 
Additionally, other UAV’s in procurement by the Department of 


the Defense (DOD) could be studied at NPS. 


pi BACKGROUND 


A. HISTORICAL DEVELOPMENT 

The history of flight testing would in itself make an 
alluring and fascinating book; this cursory summary of 
parameter estimation and flight testing reviews but a small 
fraction of the significant events in the history of flight 
testing. The majority of the historical content was found in 
the opening remarks given by Herman A. Rediess [Ref. 5] ata 
1973 symposium on parameter estimation techniques. 

One of the first test programs to obtain quantitative 
measurements of aircraft aerodynamic characteristics in flight 
was reported by Warner and Norton in 1919 [Ref. 6]. Tests 
were conducted on two Curtiss "Jenny" JN-4H biplanes at 
Langley Field, Virginia. Lift and drag coefficients were 
estimated by measuring airspeed, pitch attitude and engine 
speed in flight and assuming certain engine’ thrust 
characteristics. This specific flight test was a meager 
beginning for in-flight testing, but today it is very apparent 
that in-flight testing is a vital requirement. A 1933 report 
by Soulé and Wheatly [Ref. 7] is thought to be the first 
report to have determined and compared major longitudinal 
stability and control derivatives from flight test data with 


results acquired through theoretical predictions. The 


airplane was the single engine Doyle O-2. The analysis used 
Simplified models, solving for one parameter at a time while 
assuming values for the other parameters based upon wind 
tunnel data or other flight tests. Early in the 1950’s a 
major advancement in parameter estimation was achieved by 
Shinbrot [Ref. 8] using least squares curve fits between the 
equations of motion and flight data. A considerable drawback 
at that time was the extensive calculations required for this 
approach. These calculations, of course, were completed 
entirely by hand as the digital computer was not yet 
available as an engineering tool. 

Significant improvements were further realized in the 


later 1950’s, and throughout the 60’s and 70’s, due to: 


1. The availability of the digital computer. 


2. The progress in the technical disciplines of system 
identification and numerical analysis. 


3. The availability of high speed automatic data acquisition 
systems. 

Again, an excellent overview of the evolution of parameter 

estimation techniques up to approximately 1970 is contained in 
Reference 5. 

There are numerous methods for extracting the stability 
and control derivatives from flight-test data that have been 
developed and tested since the early 50’s. Each starts with 
equations of motion and essentially attempts to curve fit 


calculated results to the flight-test data by adjusting each 


of the derivatives or coefficients in the math model. Some, 
but certainly not all, techniques that have been used with 
success include: Ordinary Least Squares, Weighted Least 
Squares, Deterministic Least Squares, Maximum Likelihood, 
Statistical Linearized Filter, Extended Kalman Filter, 
frequency domain methods, and an older technique used in the 
1950’s, called analog matching. This last technique was a 
manual curve fitting method using an analog computer in which 
the results were very much dependent upon the skill of the 


operator. 


B. MODERN DEVELOPMENT 

Major contributions to aircraft parameter estimation since 
the mid 1960’s have been made at two NASA facilities, the NASA 
Ames Research Center’s Dryden Flight Research Facility and the 
NASA Langley Research Center. The parameter estimation 
contributions from these facilities have been made primarily 
through the work of Lawrence Taylor, Kenneth Iliff, Richard 
Maine and James Murray. 

Taylor and Iliff developed a Newton-Raphson parameter 
estimation program in 1966 based upon the theoretical work of 
Balakrishnan [Ref. 9]. The program used a modified Newton- 
Raphson algorithm to effect the maximum likelihood technique 
for estimating stability and control derivatives. This 
program underwent a gradual evolution during its application 


feet. 10]. The outcome in 1973 was a program named MMLE 


(Modified Maximum Likelihood Estimator) which used the same 
basic algorithm (Newton-Raphson) but incorporated features 
useful for processing large amounts of flight data. This 
program was widely circulated among industry and government 
agencies and as of 1980 had been used to analyze over 35 
different aircraft at the NASA Dryden Flight Research Center 
alone. [Ref. ll:p.2] 

Development of MMLE3 (Modified Maximum Likelihood 
Estimation program version 3) was completed by Maine and Iliff 
in 1980 in response to a requirement for a more versatile 
parameter estimation program. MMLE3 had two advances over 
MMLE: more flexibility in defining the equations of motion 
(although still linear); and the capability for estimation in 
the presence of state noise, also called process or input 
noise, a good example being atmospheric turbulence. ([Ref. 
ll:p. 2) Further details on MMLE3 will be addressed later. 

In 1987, a new parameter estimation program to accommodate 
nonlinear models was developed at NASA Dryden by Maine and 
Murray [{Ref. 12]. This parameter estimation program, named 
pEst, supports nonlinear models, and thus aircraft dynamic 
behavior can be tested at extreme flight conditions (high AOA) 
and for unique configurations (oblique wing, etc.).[Ref. 2] 

The basic concepts of aircraft parameter estimation 
techniques have remained unchanged for over two decades and 


are shown in Figure 1 [Ref. 5:p.14]. These concepts include: 


Noise 


Control 


Input 
Test A/C eC Measured Response 


Gontral : 
Input Aircraft 


Computed Response 
Math P P 


i Lmati Criterion Response Error 
Algorithm Function 


"Best "| Estimates 
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Figure 1. Basic Concepts of Contemporary Parameter 
Estimation Techniques 


(1) the mathematical model, (2) the data acquisition system, 
(3) the estimation algorithm, and (4) the required test 
inputs. Each of these elements will be discussed later in 
more detail. 

A Maximum Likelihood (ML) estimator was chosen as our 
parameter estimator for the following reasons. First, this 
method has become and still is "accepted as the standard 
approach to determining aircraft stability and control 
derivatives from flight data" [Ref.13:p.558]. Furthermore, 
the flight regimes of the test vehicles at NPS, at least 
initially, are expected to be well within the region where a 
linear math model will provide accurate parameter estimations. 
Furthermore, a PC compatible ML estimator program was chosen 
because of the PC’s flexibility and availability at NPS. This 


combination appeared ideal for use in analyzing the ongoing 


NPS flight research and also for use in the classroom as an 


enhancement to existing teaching aids. 
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III. MATHEMATICAL MODEL 


In this section, the linear equations of motion used to 
describe a typical manned or unmanned aircraft will be 
developed. Similar developments are shown in most aircraft 
dynamics textbooks. References used in this development were 
Airplane Flight Dynamics by Roskam [Ref. 14], classroom notes 
by NPS Professor Louis Schmidt [Ref. 15], the textbook Flight 
Stability and Automatic Control by Nelson [Ref. 16] and 
Reference 2. In creating simulated data a gust or turbulence 
model was used and that too is developed in this section. 
Lastly, a discussion of the observation equation corrections 


1s presented. 


A. MODEL EQUATION DEVELOPMENT 
This development begins with Newton’s linear and angular 


momentum equations: 


F-—S (mi) (la) 
— ad — 

wages 1b 
M ae ‘ft (1b) 


Where the force, F, is the sum of the externally applied 
forces and the moment, M, is the sum of the applied moments 
about the center of gravity (cg). The use of non-rotating, 


earth reference coordinates for equations la and 1b is 


Gat 


unwieldy for two reasons. First, required measurements are 
predominately made in the rotating body axis system; and 
secondly, but of more significance, the inertia matrix or 
tensor is a function of time in the non-rotating system. 
Therefore, the axis system chosen is the standard body axis 
system, shown in Figure 2 [Ref. 15]. The body axis system is 
used by Iliff and Maine in Reference 2, but the reader should 
be advised that the equations of motion are also at times 
derived using the stability axis system. Reference 14 has a 
good description of the differences between the two axis 
systems. The origin iS positioned at the vehicle’s cg with 
the X-direction pointing out the nose of the aircraft, the Y- 
direction out the starboard wing and the Z-direction out the 
bottom of the aircraft. Transforming equations la and 1b into 


the rotating body axis system is done below: 


Poo (mV) +x (mV) (2a) 
a= 2 (H) +6xH (2b) 
OG 


where the angular momentum, H, is the inner dot product of the 
aircraft mass moment of inertia matrix, [I,,.], with the angular 
rotation vector . The inertia matrix in the body axis 
coordinate system is not a function of time as it is in the 
earth reference. The above equations are vector equations and 


can be written into scalar components. 
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The components in the body axis for the @ vector are p, 
q, and r for the roll, pitch and yaw rates, respectively. The 
components in the body axis for V are u, v, and w for the X, 
Y, and Z components of velocity, also shown in Figure 2. 

The applied forces on the aircraft can be broken down into 


aerodynamic, gravitational and thrust components. (The thrust 
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Figure 2. Body Axis System and Notation 


1s assumed to act along the X-body axis.) 
The moments can be broken down into just aerodynamic 


components, since the thrust is assumed, and gravity forces, 


3 


by definition, act through the cg, and their moment 
contributions are zero. The moment components are then due to 
the aerodynamic forces and are shown in Figure 2 as L, M, and 
N. 

Since the gravity force components in the body axis system 
depend on the orientation of the airplane relative to earth- 
fixed coordinates (assuming a flat non-rotating earth), it is 
necessary to describe the orientation of the aircraft relative 
to the earth. Euler angles are introduced to accomplish this 
transformation between coordinate systems. The Euler angles, 
Y, ©, and M, are three consecutive rotation angles needed for 
the transformations from one axis system to the other. The 
angle YW is referred to as the yaw angle. The angle © is the 
pitch angle. The angle M is the bank or roll angle. 

The aircraft is further assumed to be rigid; that is, the 
mass particles remain at constant distances from each other. 
The X-Z plane is assumed to be a plane of symmetry and thus, 
the products of inertia I,, and I,, are zero. In this model, 
the rotating engine parts and sloshing fuel are being ignored, 
and over short periods of time when data are collected, the 
mass of the aircraft is considered to be constant. 

The force component equations that are derived with the 


above assumptions are: 
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F=m(u+qwtrv) (3a) 


FY=m(v+ru+pw) (3b) 


F=m(w+pv+qu) (3c) 


Z 


and the moment component equations are: 


aid = ea 1) ~DO1 | (4a) 
Minit och raro Tere ea)" aT (4b) 
fee ee 1) trl) (4c) 


Where the left hand sides of the above force and moment 


component equations are: 


F,=qsC,-mgsin6+Thrust (5a) 
Fy,=gsC,+mgsingcos® (5b) 
F,-asC,+*mgcos$cos6 (Se) 

M,.=qsbCc, (6a) 


eS 


M,=qsccC,, (6b) 


M.=qsbC,, (6c) 


where G 1s the dynamic pressure and q is the pitch rate. 

The vehicle is assumed to operate at small side slip 
angles and small perturbations around ae steady state 
condition. The perturbations are assumed to be small in order 
that the sines and cosines of the disturbance angles are 
approximately the angles themselves and one (1) respectively, 
and the products and squares of the products are negligible 
when compared to the quantities themselves. This 
approximation is termed small perturbation theory and permits 
the equations of motion to be decoupled and linearized into 
two smaller subsets: lLateral-Directional and Longitudinal. 

It 1s often times more convenient to have the equations in 
terms of a (angle of attack), $ (sideslip angle), and V than 
in terms of u, v, and w. These angles are usually measured 
directly vice the velocity components. In the transformation 
of the equations of motion into equations with @ and f, the 


following can be noted: 


a=tan!(—) (7a) 
ul 
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B=tan*(—) (7b) 


and if the angles are small, then 


o = (8a) 
U 

gn” (8b) 
U 
V 

oo 8 

B - (Se) 

Gees (8d) 
u 


These equations will be used in the next sections to construct 


the basic aircraft model. 


B. SIMPLIFIED LONGITUDINAL EQUATIONS 

The longitudinal set of equations pertains to rotation or 
moments about the Y axis (4b and 6b) with translation or 
forces along the X and Z axis (3a, 3c, 5a and 5c). With the 


substitution of 8a and 8b from above, and also with the use of 
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small perturbation theory, the following Simplified 


longitudinal equations are formulated: 


&= IE cig+Z (9) 
mu u 
oF (10) 
Y 
Ne (11) 


where the lift is approximately parallel to the Z-axis, so 


Gre zie. (12a) 
Ch=C, a+C,, O,4+C, (12b) 

on Cc 
Cy On aoe Cr (12c) 


Equations 12a and 12b are then substituted into equation 9 
while equation 12c is substituted into equation 10. The 
result expressed in state-space format and using dimensional 


derivatives follows: 
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iieequation 13a the M, term includes the M; term. The measured 
variables are V, a, q, 9, A, (normal acceleration in "g’s") 
and 6, during the maneuvers. The magnitude of the velocity, 
V, is approximately equal to u at small QQ and has been 
substituted for u in the above state space equation. The 
dimensional derivatives formulation is shown in the list of 


symbols. 


C. SIMPLIFIED LATERAL-DIRECTIONAL EQUATIONS 

The lateral-directional equation set pertains to the 
rotation about the X and Z axes (4a, 4c, 6a and 6c) and 
translation along the Y axis (3b and 5b). These equations are 
used to derive the lateral-directional state space 
representation. Small perturbation theory, and the use of 
equations 8c and 8d, are used to obtain the following lateral- 


directional equations: 
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Substituting equations 14a, 15a, and 16a into 14b, 15b and 16b 
respectively, and again expréssing these equations along with 
equation 17 =.in state-space format using dimensional 


derivatives the following is obtained: 
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The measured variable are the 6, and 6, control deflections, 
B, p, x, ®, A, and v. 

The two state space representations (longitudinal and 
lateral-directional) form the mathematical model used with the 
parameter estimation program to estimate the stability and 


control derivatives. 


D. TURBULENCE MODEL 

In preparation for using actual time history data, 
Simulated longitudinal and lateral-directional data were 
created. In creating the simulated data it was desired to 
match the expected phenomena or real world effects that would 
be encountered. Thus turbulence (state noise) and measurement 


noise were selectively added to the simulated data by the 


od 


user. This section describes the model used to generate the 
turbulence or state noise. 

The development begins with the application of the Dryden 
gust model [Refs. 15 and 18]. The turbulent effect can be 
added to both the longitudinal and lateral-directional 
equations. The development for the longitudinal case will be 
shown below and can be extrapolated in a straightforward 
manner for the lateral-directional case. 

In the Laplace domain the vertical gust velocity transfer 


function is 








Wy (8) =VE-S* (5) (19) 
where 
b= = (20a) 
Lag (20b) 
. i (20c) 


Ze 


and L, is the scale of the turbulence, 0, is the rms value of 
the turbulence and TN] 1S a zero mean white noise input. The 
values used for L, and 0, equated to a turbulence level 
between light and moderate, and can be adjusted if desired. 
Equation 21 was obtained by transforming equation 19 from 


the Laplace domain into the time domain and dividing by u: 


Wo = +2h6,+h?04= VX (by +H) (21) 


Meera, 1s the & perturbation) attributable to the gust. 


Equation 21 is converted into the state space format by 


letting 
Ty Ms (22a) 
= WES 
2,78, aie (22b) 
Z,27,+VKy (22c) 
u 
#,=-242,-072,+V Ey (b-2d) (22d) 


me Lollows that 
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This result can be combined with the longitudinal state space 
equation 13 (Lateral-Directional equation 18) to yield 


equation 24a: 
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Similarly, the state space equations for the lateral- 
directional model with turbulent noise can be developed and 


are shown below: 
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E. MEASUREMENT CORRECTIONS 

The observed or measured data must be corrected for 
measurement errors caused by the positioning of the sensors. 
The aircraft equations of motion were developed for the 
masecraft cg. Therefore, measurements taken by sensors not 
physically located at the cg require corrections to reference 


the values to the cg location. This subsection will document 
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the correction equations used in this analysis. The simulated 
data need not be corrected since these data were manufactured 
at the cg of the aircraft. However, in anticipation of actual 
aircraft flight test parameter estimation, sensor position 
corrections were implemented into the programs used to analyze 
genuine flight data. 
1. Longitudinal 
The longitudinal a and A, data are capable of 
corrections for sensor position displacement from the cg. The 


Q& was corrected for the X-coordinate probe position forward 


(+) or aft (—-) of the cg (Xf) as shown below: 
X 
Og" 8probet > F (26a) 


A correction for the upwash angle @, at the probe was not 
taken into account as it 1s assumed to be small or previously 
accounted for when the sensor was calibrated; a more complete 
discussion on corrections 1s contained in Reference 2. The 
normal acceleration, A,, was corrected for both X,, (fwd +) and 


Zan (Gown +) displacement from the aircraft cg as shown below: 


Xx Z 
A =A a ee (26b) 
g g 
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2. Lateral-Directional 


The sideslip angle, $, and the lateral acceleration, 
A,, were corrected for their sensor positions and the 


correction equations are shown below: 


Xp 
B.g*B prope” —> I (27a) 


and 


A, =A, Sata (27b) 
cy acc g g 


foemevalues for x, 4%, and 2, are again defined positive 


forward and down analogous to the body coordinates. 
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IV. MAXIMUM LIKELIHOOD PARAMETER ESTIMATION 


A. THEORY 

The concept of parameter estimation, as discussed earlier, 
can be defined quite simply in general terms. The system, in 
this case a UAV or some other testbed whose parameters are to 
be estimated, is assumed to be described by a set of linear 
dynamic equations, a mathematical model, which was defined 
earlier. To determine the values for the unknown parameters 
the system is excited by an input. The input and the system's 
actual response are measured. The values of the system 
unknowns are then calculated based on the requirement that the 
model response to the input match the actual system response. 
Complications to this simplified explanation arise when the 
following are considered: 

1. Measurement noise —- perfect measurements are unattainable 
with any sensor. 


2. State noise - the aircraft or system is being excited by 
unmeasured sources such as atmospheric turbulence. 


3. Modeling errors - exactly describing the physical system 
with simple, especially linear, dynamic equations is very 
unlikely. 


In fact, 1f the above complications were not present, the 


exact values of the unknown parameters could be found and what 
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is termed parameter "identification" vice "estimation" would 
be our accomplishment. 

The common approach for handling the modeling error is to 
ignore it and let the error be treated as measurement or state 
noise or both. I1iff and Maine state that "this procedure is 
not rigorously justifiable, but combined with a carefully 
chosen model, it is probably the best approach available." 
feet, 19:p. 2) The Modified Maximum Likelihood Estimation 
algorithm version 3 (MMLE3) program was structured to take 
into account the presence of state and meaSurement noise. 
The information in the following section is a compilation of 
information on MMLE parameter estimation from References 2 and 


ieenrough 22. 


B. MODIFIED MAXIMUM LIKELIHOOD ESTIMATION 

In this section the Modified Maximum Likelihood Estimation 
algorithm version 3 (MMLE3) is presented. The first step for 
parameter estimation then, as mentioned above, is to model the 
system accurately. That model was developed in the previous 
chapter. The aircraft equations of motion define the system 


model and can be expressed in the following state space form: 


X(t) =X (28) 


RiGee wat) + BujiGe) +hw(t) +Dy (29a) 
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y(t) =Cx(t) +Du(t) +b, (29b) 


Z(k) =y (k) +Gv(k) (29) 


where x(t) is the state vector (x, being the initial state), 
u(t) is the control input vector, and y(t) 1s the prediction 
or model output vector. Matrices A, B, C, and D contain the 
unknown system parameters, which in this case are the 
stability and control derivatives. Matrices F and G represent 
the covariance matrices of the state and measurement noise 
respectively. The measured response vector, z(k), is sampled 


at N discrete time points (k=1,...,N). 


The state or process noise, w(t), 1s assumed to be zero- 
mean, white Gaussian and stationary. The measurement noise, 
v(t), 18 assumed to be zero-mean Gaussian noise with identity 
covariance. 


The complete unknown parameter vector to be estimated is 


then given by: 


€=(H 7A?) Deby) (30) 


where H represents the unknown parameters in the matrices A, 
B, C, and D; 2X represents the unknown elements of the F 
matrix; and b, and b, represent the unknown biases of the state 
and model output equations respectively. The [-]’ indicates 


the transpose of a matrix. 
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The maximum likelihood estimates are obtained by 
minimizing the negative logarithm of the likelihood function. 
The likelihood function, J, 1s a function of the difference 
between the measured and computed time histories. The 


likelihood function is: 


N 
= . Tp-1 z N 31 
WGra) 2 dy |2 z(k) -y(k))7R+[2(k) y(k)]+>1n|R| (31) 


where Ris the innovation covariance matrix. The innovations 
or residuals are [z(k)-y(k)]. In order to obtain the 
predicted state variables it 18S necessary to use a state 
estimator. The Kalman filter, which is an optimal linear 
state estimator, is used for this purpose. The Kalman filter 
consists of a prediction step and a correction step [Ref. 


20:p. 12] for equations 29a and 29b and is shown below: 


X(k+1) =>2(k) +WBu(k) +b, (32) 
y (k) =C&(k) +Du(k) +b, (33) 
R(k) =¥(k) +K (Zz (k) -y(k) ] (34) 
u(k) = = + (u(k+1) +u(k)] (35) 


sal 


where ~ (tilde) and ~ (circumflex) denote the predicted and 
corrected state variables respectively. K represents the 
Kalman filter gain matrix. The state transition matrix is ®M 
and its integral is Y. 

1. Cost Function Minimization 

The maximum likelihood estimates for the unknowns are 
found by minimizing the negative logarithm of the likelihood 
function, J(€,R). The negative logarithm of the likelihood 
function is often called or referred to as the cost function. 
This minimization is done by using a modified Newton-Raphson 
technique which iterates on the vector of unknowns, ¢, with 
each iteration providing a new estimate of the unknown vector. 
These new estimates update the math model coefficients, 
providing a new calculated response and a new response error. 
This iteration process is continued until the convergence 
criterion is satisfied. 

"The maximum likelihood estimation method has the 
desirable characteristics of yielding asymptotically unbiased, 
consistent and efficient estimates [Ref. 19:p.3]." 

2. A-Priori Weighting 

The MMLE3 algorithm allows for the use of a weighting 
function to account for prior ’engineering’ knowledge of the 
aircraft parameters. This prior knowledge can be obtained 
from other test cases, wind tunnel measurements, or indepth 


analytic analysis. A table of the relative importance and the 
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prediction accuracy of stability derivatives using theoretical 
methods is contained in Reference 14, page 236 and can be used 
as a guide, if desired, to weighting the initial parameter 
estimates. 

The a-priori information can assist the algorithm in 
converging, but caution should be used, as the weightings can 
prejudice the answers toward the analyst’s own values 
feer.22:p.ST-8]. 

3. Estimate Uncertainty 

The use of the Cramér-Rao bounds with the maximum 
likelihood estimator can also provide a measure of the 
relative accuracy of the estimates. Each parameter bound 
gives an approximation to the standard deviation of the 
estimates. It 1s important to recognize that these bounds 
are, in fact, the lower limits for the standard deviation, 
meaning that the standard deviation value is at least as large 
as the Cramér-Rao bounds [Ref. 21:p.12]. More details on the 


Cramér-Rao bound is contained in References 2, 20, 21 and 22. 
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The State Space Identification Toolbox (SSIT) for the 386- 
MATLAB personal computer program implements the MMLE3 
algorithm. MATLAB is a registered trademark for matrix 
oriented software distributed under license agreement by The 
Mathworks, Inc. Reference 21 18 a report of the results of a 
study comparing the mainframe based MMLE3 program and the 
MATLAB SSIT implementation of MMLE3. The analysis indicated 
that the PC version results were "generally well within the 
uncertainty levels of the mainframe parameter estimates [Ref. 
Zi. es The use of MATLAB Software will be discussed in 


the following chapter. 
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Vv. APPLICATION 


A. DATA ACQUISITION 

The data acquisition system is an important part of 
dynamic stability and control testing. The more information 
known on the details of the entire data acquisition system, 
the greater the probability that the test results can be more 
precise. With few details known about the data, often times 
only gross characteristics of the aircraft can be determined. 
The details essential for a complete analysis of the data 
should include how the data were filtered, digitized, time 
tagged, transmitted and recorded. The complete analysis of 
the data acquisition system should start at the beginning with 
the sensors and continue through to the final recorded data 
Pproauct. 

Sensor calibration errors, temperature effects, added 
noise from aircraft vibration, recorders and transmitters are 
but a small portion of the circumstances that should be 
reviewed. Common recording systems, their advantages and 
disadvantages, plus other issues relevant to the entire data 
acquisition process are discussed in greater detail in 


Chapters VII and VIII of Reference 2. 
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B. INPUT SELECTION 

The selection of the control inputs for use in parameter 
estimation must take into account the pilot’s acceptability 
and safety of flight concerns as well as the model validity 
considerations. 

References 1, 2, 22 and 23 detail various methods for the 
input design employed in aircraft parameter estimation. 

One specific requirement is that the controls applicable 
to the specific model need to be exercised, such that the 
aircraft modes are excited. Control inputs which are near the 
frequency of the excited mode usually provide the best 
results. This is because at these modal input frequencies, 
the largest aircraft response for a given input usually occurs 
and provide the estimator significant data as compared to the 
noise (gust and measurement noise). Judgment must be used 
when selecting the control inputs to insure flight safety and 
to avoid responses that exceed any preset magnitude 
restrictions. In the assumed linear model, for example, as 
the response magnitudes exceeded the small perturbation 
assumptions (10-15 degrees), the final parameter estimation 
results can be expected to worsen. Likewise, the signal-to- 
noise ratio of the data improves proportionally with the 
magnitude of the response; thus there is a trade-off in the 
development of the aircraft control inputs. Reference 23 


discusses a method to optimize the control inputs while 
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accounting for specific restrictions or trade-offs as 
mentioned above. 

The control inputs used to generate the simulated data for 
this study were elevator, rudder and aileron ramped doublets. 
This type of control input was selected for two reasons. 
First, this input is truly representative of actual pilot 
inputs (Impulses and step inputs are not physically 
realizable), and second, it can be easily adjusted in 
Magnitude and frequency as necessary for different aircraft. 
The previously mentioned references provide additional 


information on the specifics of designing control inputs. 


C. MATLAB 

MATLAB is a commercially available software package for 
scientific and engineering applications. The program 
integrates numerical analysis, matrix computation and graphics 
into a relatively simple environment without the need for 
traditional programming knowledge. MATLAB has specialized 
toolboxes for added capabilities. In this study the Control 
Systems Toolbox and the State-Space Identification Toolbox 
(SSIT) in addition to the main 386-MATLAB program were used. 
The SSIT implements the MMLE3 algorithm, discussed previously 
in Chapter IV. 

Again, the use of the MATLAB based program is desirable at 
NPS because it operates in a familiar PC environment, the 


importation of data is relatively straightforward and easily 
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accomplished, knowledge of a formal programming language is 
not required, and the plotting and hard copy functions are 
easy to use and manipulate. Furthermore, during the required 
basic dynamics and linear systems courses, students at NPS 
use MATLAB extensively as an instructional and problem-solving 
aid. 

1. M-files 

MATLAB is capable of executing sequences of commands 
stored in files, called M-files or macros, from a single-line 
command. The M-files have a file type of .m and consist of a 
sequence of normal MATLAB statements that can include the 
execution of other M-files. Major benefits of the M-files 
are: repetitive or long sequences of commands can be 
automated; new functions can be created by the user for a 
specific need; and the .m files are ASCII type format and 
easily edited. 

The following sections will describe the M-files 
created during this study for use in implementing the PC 
MATLAB parameter estimation program. It is noteworthy that 
the SSIT is itself an M-file (mmle.m). More detailed 
information concerning MATLAB is contained in the MATLAB 
user’s guides References 22 and 24. 

a. Simulated Data 

Simulated time history data were created using 


MATLAB M-files. This simulation was done in preparation for 
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using actual aircraft flight data in the parameter estimation 
program. The aircraft models (longitudinal and lateral- 
directional) used to create the data were previously 
discussed, as was the turbulence model. The M-files described 
in the following sections are contained in Appendix A. 

(1) Longitudinal 

The simulated longitudinal data were created 
using the M-file LONGDAT.M (see Appendix A). This M-file is 
extensively commented which allows the user to understand the 
program and change or adjust certain parameters as necessary, 
such as turbulence level, measurement noise level, the 
elevator input amplitude, and period, or to design a 
completely new input. 

Aircraft derivatives and other physical data 
are necessary to create the simulated data. These data are 
stored into MATLAB data files for a small number of specific 
aircraft. These data files are .mat type files. The user 
can select one of these aircraft or input the data for an 
aircraft of his or her choosing. If a new aircraft is 
selected, the data required by the program are interactively 
requested using input commands. The data are then stored and 
available ina .mat data file for later use. The storing of 
these data into accessible files eliminates the need to 
reenter the data every time additional simulated data are 


desired for the same type aircraft. Alreraft data from 
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Reference 16 were initially entered for the following 

aircraft: NAVION, A4D, F104A, JETSTAR, and B747. The Pioneer 

UAV data were also entered and were obtained from Reference 4. 

The general input requirements for the LONGDAT.M macro are: 
1. The aircraft physical data, if not uSing an aircraft with 
previously saved data 


2. The selection of either adding or not adding state and 
measurement noise to the data 


3. The flight specifics: velocity, pressure altitude, and 
outside air temperature 
The M-file uses a ramped elevator doublet for 
the input to the aircraft math model, which then produces the 
output. The Simulation can be done either with or without 
noise. When noise is selected, in addition to the state 
turbulence noise added, a uniform measurement noise is also 
added to the outputs. The a noise added to the 
outputs @ and © is zero mean, +% degree maximum value, while 
the measurement noise added to the output gq 1s zero mean, +2.5 
degrees per second maximum value and the measurement noise 
added to the normal acceleration, A,, 18 zero mean, +.01 g 
maximum value. 
The final output time histories for the user 
include 6, a, q, © and A, plots displayed on the PC monitor, 
with data files saved containing the simulated time history 


data. The data files that are created by LONGDAT.M are then 
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available for the parameter estimation M-files, which will be 
discussed later. 
(2) Lateral-Directional 

The simulated lateral-directional data were 
created using the macro LATDIR.M. This file is similar to 
LONGDAT.M but uses the lateral-directional math model. The 
model has two inputs, 6, and 6,. The necessary constants and 
Stability and control derivatives again are stored into a .mat 
file as was done for the longitudinal case. The noise is 
essentially the same as in the longitudinal model except that 
the turbulent gusts are caused by a perturbation in the side 
slip angle, B,. The uniform measurement noise that is added 
to the output angles (f and ®) is zero mean, +% degree maximum 
value, while the noise added to the output angular rates (p 
and r) is zero mean, +2.5 degrees per second maximum value, 
and the noise added to the lateral acceleration is the same as 
in the longitudinal case. 

The outputs are time history plots and data 
files of rudder and aileron inputs, 6, and 6,; side slip, 6; 
roll rate, p; yaw rate, r; roll angle, @®; and lateral 
acceleration, A,. The data files that are created by LATDIR.M 
are now available for the parameter estimation M-files, which 


will be discussed next. 
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b. Parameter Estimation 

Once the required time history data are available, 
either created by simulation or acquired from actual flight 
tests, the parameter estimation algorithm is used to calculate 
the ‘best’ estimate of the stability and control derivatives. 
The process of executing the MATLAB SSIT parameter estimation 
program was accomplished with four basic M-files. The four M- 
files were developed so the execution of the SSIT M-file 
(mmle.m) would appear transparent to the user. Modifications 
to the basic four M-files were necessary for each of the 


following cases: 


1. Simulated longitudinal data 
2. Simulated lateral-directional data 
3. Actual longitudinal flight data 


4. Actual lateral-directional flight data 


These M-files are included in Appendix A. The four basic M- 
files used with the simulated data will be discussed followed 
by the changes needed for the actual flight data cases. 

The four basic M-files used in each of the above 
four cases were designed to simplify the execution of the SSIT 
M-file (mmle.m). This M-file arrangement is shown by the 


block diagram in Figure 3 below: 
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USER 
INTERFACE 


MATLAB 
NPSMMLE SSIT 


(MMLE.M) 


NEsicsS 


NPSINIT MLEPLOT 


OUTPUT 





Figure 3 Block Diagram of M-file Arrangement 


The four cases mentioned above each have the four 
basic M-files shown in the above block diagram. The above 
case number is included in the M-file name used for that case 
to distinguish the files between the different files. Four 
files for each of the four cases equate to 16 different M- 
files. Each set of files is a variation of the basic four M- 
files. 

The only M-file the user initiates is the macro 
Gemed NESMMLE (the 1s where the case number from above is 
included with the M-file name) shown in Figure 3. The other 
macros are "called" in a manner Similar to that for a 
subroutine call in FORTRAN programming. 

NPSMMLE is initiated by the user and it acts as 
the link between the SSIT, the other macros and the user. The 


one exception is that the time between data points, dt, must 
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be edited in NPSP2SS_.M for cases 3 and 4. The functions of 
NPSMMLE include the following: 
1. It initiates the MATLAB diary function which saves a copy 
of the output (no plots) for subsequent printing 
2. It "calls" the NPSINIT macro (discussed later) 


3. It arranges the time history data into the required 
format 


4. It establishes and makes known the function file to the 
SSIT, which converts the parameter vector into the user 
defined state space description of the model 


5. It defines additional information to the SSIT which is 
discussed in detail in References 21 and 22 


6. It "calls" the SSIT parameter estimation program (mmle.m) 


7. It formats and displays the final numerical results to 
the monitor 


8. It "calls" the MLEPLOT_ macro which graphically displays 
the results 
The NPSINIT_ macro is "called" by the NPSMMEEE 
macro and performs the following functions: 
1. It loads the time history data file supplied by the user. 
This data file is either obtained by simulation or from 
actual flights. 
2. It requests the time interval between data points. 
3. It requests the vehicle’s physical attributes, initial 
parameter estimates and the flight conditions. 
All of the above information is saved and then made available 
for the SSIT when it is reloaded into the MATLAB working 


environment by NPSMMLE . 
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The NPSP2SS_ macro is a function file used by SSIT 
during the parameter estimation process. The file converts 
the parameter vector into the user defined stated space 
description of the system. Details of the P2SS M-file are 
contained in References 21 and 22. 

After the SSIT has been "called" and successfully 
run, the numerical results are output to the analyst by way of 
foeeeC monitor. The MLEPLOT macro is then initiated and its 
function is to display the results graphically to the user. 
These plots are shown on the monitor and saved as MATLAB .MET 
files. These .MET files are high resolution graphics files 
that may be used later for printing graphics hard copies. 

The macro file versions numbered 1 or 3 are for 


use with the longitudinal data, simulated and actual data 


respectively. The versions numbered 2 or 4 are for use with 
the lateral-directional data, Simulated and actual 
respectively. 


The differences between the versions 1 and 3 and 
2 and 4 are relatively small. The biggest differences are 
that when the actual data (3 and 4) are used, the user is 
asked for sensor position corrections and queried whether an 
a-priori weighing vector is to be used during the parameter 
estimation process by the SSIT. If selected, the weighting 


input is a vector of the variances for each of the parameters 
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to be estimated. Therefore, the higher the number, the less 
the weighting afforded that parameter and vice versa. These 


M-files are included in Appendix A. 
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ne RESULTS AND DISCUSSION 


A. SIMULATED DATA 
1. Application 

The MATLAB M-files LONGDAT.M and LATDIR.M were used to 
develop the required simulated data representing time history 
responses. These data included both state and measurement 
noise. Qualitatively, these data, plotted in Appendix B and 
discussed below, compare favorably with actual time history 
data shown in many references. Thus, the model chosen for the 
Simulation was assumed as an adequate mathematical 
representation of the physical system within the region of 
applicability, this region being the area or flight regime 
satisfying the restrictions placed upon the math model during 
development. Additionally, with the simulation, there is the 
capability for the user to modify the control inputs and noise 
magnitudes thereby enhancing the use of these programs as an 
iistructional aid. Various aircraft types with differing 
inputs and atmospheric conditions potentially could be 
examined. 

The parameter estimation of the simulated data was 
quite accurate with some exceptions. These exceptions will be 
discussed later. The accurate parameter estimates, for the 


given model, validate the MMLE methodology in the presence of 
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both state and measurement noise. The results show the 
relative insensitivity of the MMLE method (implemented by the 
MATLAB SSIT) to noise when the modeling and the excitations 
(inputs) are chosen with care. 

Longitudinal and lateral-directional parameter 
estimation examples, plots and quantitative output using the 
Simulated data are shown in Appendix B. These simulated 
examples are for the A-4D and Navion aircraft and the PIONEER 
UAV. The SSIT quantitative results are tabulated with the 
initial input parameter estimates and the "truth" Or 
underlying parameters (parameters used to generate the data) 
for comparison. The tabulated SSIT quantitative results are 
presented in column format with the column headings as 
follows: pid, parameter id number; p(pid), final parameter 
estimate; pref, initial input reference parameter; cramer, 
Cramer-Rao bounds; 2fcramer, two times a corrected, filtered 
Cramer-Rao bound; insens, insensitivity, the change in that 
parameter required to move from the maximum likelihood value 
to the edge of the confidence ellipsoid. For a single 
parameter model the insens value is the Cramer-Rao bound. 

2. Longitudinal 


a. A-4D 


The SSIT parameter results for the A-4D, due to 


the elevator control input shown in Figure B-1l, are presented 
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Pier igures BeZ wazehrough B5B-5. These are the input and 
subsequent response of the A-4D. The a, gq, 9, and 
acceleration responses correlate well with the truth data used 
to generate the simulated data. All estimated parameter 
values for Cyg,r Char Cugr Crsee, and Cys, Shown in Appendix B are 
within the 2fcramer bound. This bound represents the 95 
percent probability ellipsoid for the parameters. Close 
inspection of the estimates show that the greatest deviations 
eee tor C,, and C, and both are within four percent of the 
actual underlying values used to create the simulation. 
b. Navion 

The SSIT results for the Navion are shown in the 
response plots, Figure B-7 through B-10, due to the elevator 
control doublet input shown in Figure B-6. Again, as was seen 
in the previous A-4D results, the correlation of the estimated 
aircraft response plots with the underlying truth derivative 
feeeenses are good. The largest errors are in the C,, and C,,, 
derivatives, and they are 12.4 and 7 percent respectively. An 
Beeurate prediction of Cy using theoretical methods is 
difficult to achieve; Roskam [Ref. 14] notes an acceptable 
estimated prediction accuracy for this derivative of 20 
percent. 

Cc: UAV 
The results for the A-4D and Navion are more 


accurate than the UAV results which are shown in Figures B-11 
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through B-15. The UAV results demonstrate the importance of 
designing adequate inputs for each specific vehicle. The 
input must sufficiently excite the vehicle’s dynamic response. 
Identical elevator input doublets were used for all three 
vehicles; the UAV elevator input is shown in Figure B-ll. The 
input period is 2 seconds with a maximum amplitude of 
approximately 4 degrees. The small UAV response due to the 
small elevator control power (C,;.) might be overlooked if the 
scaling on the response plots were not closely observed. The 
responses are less than a half to a third those of the A-4D 
and Navion. These small responses equate to a lower signal- 
to-noise ratio for the parameter estimator. The response due 
to the elevator input was not much more significant than the 
response perceived by the estimator from the state noise (the 
gust) in combination with the presence of the measurement 
noise. Thus, the parameters with the exception of C, (2.4 
percent) were all in error greater than 30 percent, and C,, was 
76 percent in error from its underlying truth value. Caution 
must be exercised in choosing a proper excitation tailored for 


the particular vehicle’s response. 


3. Lateral-Directional 
The lateral-directional parameter estimates provided 
by the SSIT for all three simulated aircraft are also shown in 


Appendix B. 
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The A-4D lateral-directional inputs are shown in 
Figure B-16 and consist of a rudder doublet followed by an 
aileron doublet. The aircraft responses of 6, r, ®, p, and 
lateral acceleration due to the aileron and rudder inputs are 
shown in Figures B-17 through B-21. As can be seen in these 
plots, the correlation between the measured and estimated 
responses is good and the parameter estimates which are also 
shown in Appendix B were accurately estimated. Of the 12 
parameters all but three were inside the estimator 2fcramer 
bound. The three derivatives (i.e. G,, Cy,3, and Cys,), although 
not inside the 2fcramer bound, were very close to the bound, 
and are within 18, three and eight percent of the truth values 
respectively. Roskam [Ref. 14 p.236] indicates that the 
theoretical accuracy for values in estimating CC, is 
approximately 25 percent using analytical methods. It is felt 
that improvements in the estimation of these derivatives could 
be achieved by investigating different inputs. 
b. Navion 
The Navion results are also shown in Appendix B, 
Figures B-23 through B-27 are the response plots due to the 
rudder and aileron inputs as described by Figure B-22. The 
inputs are identical to the rudder and aileron inputs used for 


the A-4D. The Navion aircraft responses 6, r, ®, p, and 
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lateral acceleration to the inputs (Figures B-23 through B-27) 
also show good correlation and indicate that the parameter 
estimates were estimated accurately. A vertical shift in the 
estimated plot from the true response sometimes occurs and can 
be misleading. An example of this vertical shift is shown in 
Figure B-26 by the Navion roll response. This misleading 
vertical shift is caused by noise on the first data point. 
The estimated plots were originated by MATLAB at the first 
data point and any noise in that point causes a vertical shift 
of the estimated curve. A skilled analyst can adjust the 
first point to the known initial flight condition and 
eliminate the misleading vertical offset. The numeric values 
for the Navion parameter estimates are also shown in Appendix 
B and were accurately estimated. Of the 12 parameters all 
were inside the 2fcramer bound. 
C: OAV 

The parameter estimation results for the UAV due 
to the rudder and aileron doublets were also quite accurate. 
The aircraft inputs and response plots are shown in Figures B- 
28 through B-33. The identical rudder and aileron inputs used 
for the A-4D and Navion were used for the UAV and are shown in 
Figure B-28. The misleading vertical shift 1s more prominent 
in the UAV responses, especially in Figures B-30 and B-33, the 


yaw rate and lateral acceleration plots. Again, the responses 
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correlate well (except for the vertical shift) and the 
estimated parameters are all within the SSIT 2fcramer bounds. 
4. Problems 

Significant problems were encountered when trying to 
employ the SSIT on a very lightly damped or divergent system. 
This case was experienced with the F104A aircraft. buring 
these conditions the SSIT mmle.m program would halt prior to 
completion due to an error. The error displayed to the 
analyst was always a matrix singularity problem encountered 
during the SSIT computations. 

Another rarely encountered problem was the case where 
the calculated estimates were significantly in error from the 
actual parameters. This case was very obvious to the user, 
especially when the estimated plot was compared with the 
measured data. These significant parameter errors (1 to 2 
orders of magnitude) caused the estimated plot to rapidly 
diverge from the measured data. It is believed that the 
estimation program was converging upon a local minimum instead 


of the global minimum of the cost function. 


B. ACTUAL FLIGHT DATA 
1. Application 
The two actual flight tests analyzed were for the 
longitudinal cases of the F-14A and the T-37 aircraft. The F- 


14A data were acquired from the Naval Air Test Center and the 


53 


T-37 data from NASA-Ames Dryden Flight Research Center. The 
quantitative results and plots are shown in Appendix B. 
a. F-14A 

The F-14A elevator input is shown in Figure B-34 
with the a, g, © and acceleration responses shown in Figures 
B-35 through B-38. The estimated response period appears to 
correlate well but the magnitude of the estimated responses 
are less than the actual measured responses in all cases. 
This result could be due to the accuracy of the aircraft 
physical characteristics used in the aircraft model or to data 
sensor location inaccuracies. These possible problems and the 
predicted derivatives will be discussed later. The estimation 
does, however, give a general representation of the aircraft 


longitudinal characteristics. 


bi. ates 7 

The T-37 elevator input is shown in Figure B-39. 
The aircraft response is shown in Figures B-40 through B-43. 
The response plots are very accurate with the exception of a 
disparity in the © response (Figure B-42) after approximately 
the 5 second point. It is not known what caused this 
perturbation since the q (pitch rate) response correlates well 
and is the time derivative of © response. A possible cause 
could be attributed to a data acquisition problem, and quite 


possibly the cause will never be known. 
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2. Problems 

The problems encountered using actual flight data 
involved obtaining accurate aircraft physical characteristics 
and arranging the data into the desired format. Since the 
flight tests were not specifically designed and conducted for 
this study, much of the required physical data were not 
readily available and were estimated. The use of estimates 
was especially true for the F-14A data. Additionally, with no 
truth data and so many variables on the F-14A such as wing 
sweep angle, independent left and right horizontal control 
surfaces and unknown external loadings, the accuracy of the 
model used and the results are at best an approximation. 

In both cases the output parameter estimates were 
obtained without the use of the program’s weighting 
capability. The parameter estimates appear to be reasonable 
with the exception of the C,, values, which are inclined 
towards the higher side of discretion. However, with only one 
set of data (one flight maneuver) for each airplane analyzed, 
no Significant numerical conclusions can be justifiably 
deducted. Continued experience with actual flight data will 
provide a database upon which decisions of proper weighting 


values can be determined. 


5° 


A. 


B. 


VEL. CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 


The following conclusions were deduced from this study: 


1. The simulated time history data generation and its use in 
the parameter estimation program worked satisfactorily. The 
accurate parameter estimates validated the MMLE method as 
implemented by MATLAB and its relative insensitivity to 
state and measurement noise if the model and inputs were 
carefully selected. 


2. It is thought that significant benefits can be achieved 
with the simulated portion of the study as a classroom 
instructional aid in the Flight Test and Flight Dynamics 
courses at NPS. 


3. The actual flight test data appeared to produce 
acceptable results. However, many details were not known 
concerning the data and aircraft characteristics. These 
unknown details concerned such items as data filtering, 
sensor position, accurate aircraft physical characteristics, 
and time lags. Direct involvement with the flight tests, 
although not essential, would have significant benefits. 

These benefits would include knowledge of the above missing 
details, easier data acquisition into the desired data 
formats and hopefully amore in-depth and accurate analysis. 


Recommendations 


The following recommendations are offered based upon the 


above conclusions: 


1. Incorporate this study into the Flight Test and Flight 
Dynamics courses at NPS as instructional demonstrations. 


2. Continue to develop an on-site data reduction capability 
to enhance the flight test research being conducted with 
UAV’s at NPS. This development will enable the school to 
have the capability of comparing computational computer 
studies and wind-tunnel studies with flight-test results. 


56 


3. Investigate the possibility of incorporating the pEst 
non-linear parameter estimation program (Reference 12) into 
the flight test research program to cope with any future 
requirements. These applications could include helicopter 
dynamics, high-a flight regimes or a host of other areas not 
particularly well suited to linear modeling. 


4. Investigate developing a neural network parameter 
estimation capability. This capability could be used for 
development of a real-time reconfigurable flight control 
system to improve aircraft survivability. These development 
areas tie parameter estimation into two additional research 
disciplines of interest at NPS, Neural Networks and Aircraft 
Combat Survivability. 


5. Investigate assimilating the GAT-1 training device and 
the results from parameter estimation tests into a 
reconfigurable simulator for instructional purposes in 
advanced Flight Dynamics, Control and Avionics courses. 


6. Investigate the feasibility of using the MATLAB SSIT to 
accurately determine ship dynamics and perhaps find better 
ways of controlling unwanted ship motion. The improvement 
in ship dynamics could lead to improved helicopter landing 
conditions during rough sea states. 


7. Finally, continue the NPS flight test research program on 
Department of Defense UAV’s such as the Pioneer and Exdrone, 
not only to assist in their evaluation but to stimulate the 
Students’ interest with relevant and available military 
research topics. 
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APPENDIX A -- M-FILES 
A. LONGDAT.M 


clear; 

% MACRO FILE NAME >=S=S===== LONGDAT .M =======< 

% Date 1 Feb 92 

a ee eee ee ee eee ee ee SS SSS See SS OS SS ee Ie 
disp (ans) 

disp(’ ‘) 

arse MACRO TO GENERATE SIMULATED LONG. DATA a 

disp (‘ USING ELEVATOR DOUBLET WITH OR W/O NOISE i) 

disp (ans) 

disp >) 

§ewn nn ee ee eee GET AIRCRAFT TYPE TO BE USED crrrr------- 
dismas AIRCRAFT TYPES AVAILABLE oo 
disper 

disp (’NAVION A4D FI1O4A JETSTAR B747 UAV OTHER a 
di Spr a) 

disp (’ SELECT OTHER TO INPUT DATA FOR USER DEFINED AIRCRAFT =) 
diusta(2 ae 

disp (’NOTE ========> PROGRAM IS CASE SENSITIVE <========' ) 
disp (es; 

typac=input (’ TYPE IN DESIRED A/C FROM THE ABOVE LIST. Se 
Go wer nr nn nr nnn nn eee DETERMINE IF NOISE IS WANTED 
disp(’ ‘’) 

sysn=input (’ INPUT A 1 TO INCLUDE NOISE AND A ZERO FOR NO NOISE. '‘); 
ndp=201; % NUMBER OF DATA POINTS -10 SEC 

dt=.05;% TIME STEP FOR THE DATA 

amp=.07; % AMPLITUDE OF DOUBLET (RADIANS) 

period=1;% PERIOD OF DOUBLET IN SEC= PERIOD + 1 SEC 

t=[0: ndp=1)) de, % TIME VECTOR 

Simdata=zeros (ndp, 5); *% SETUP DATA MATRIX ALL ZEROS 

aa a ee eee GENERATE THE INPUT DOUBLET 


dslope=(4*amp*dt) /period; 

d1l=(0:-1*dslope:-1l*amp) ;dla=-amp*ones (1:10); 

d2= (-1*amp+dslope:dslope: amp) ;d2a=-1*dla; 

d3= (amp-dslope:-1*dslope:0); 

simdata(:,1)=[dl dla d2 d2a d3 zeros(1,ndp-(period/dt) -21)]’; 


%----------------- GENERATE THE STABILITY AND CONTROL MATRICIES 
disp(’ ‘) 
vtrue=input (’ INPUT AIRCRAFT TRUE VELOCITY ft/sec Ve 
altft=input (‘ INPUT AIRCRAFT ALTITUDE IN ite ) 
oat=input (’ INPUT THE OUTSIDE AIR TEMPERATURE Deg F ‘ys 
if (strcmp (’ OTHER’ ,typac)>0); 
typac=input (’ INPUT THE A/C TYPE ....< 6 characters | Sa 
sref=input (‘INPUT THE AIRCRAFT REFERENCE AREA £ft%*2 “ye 
gw=input (’ INPUT THE AIRCRAFT GROSS WEIGHT lbs ove 
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iyy=input (‘AIRCRAFT Iyy MOMENT OF INERTIA slug-ft%*2 ae; 
oF 


CLde=input (‘INPUT DERIVATIVE CL de 1/RAD 
CMde=input (’ INPUT DERIVATIVE CM de 1 /RAD 
ans=[{’save ‘’,typac,’ .mat;’]J; 
eval (ans) ; % SAVE THE NEW A/C DATA IN A .MAT FILE 
truth=[{CLa CMa CMq CLde CMde]; 
else 
ans={' load ’,typac,’” .mat;"]; 
eval (ans) ; 
truth=[{CLa CMa CMq CLde CMde]; 
end 
% CALCULATE DENSITY, DYNAMIC PRESSURE, AND CONSTANTS 
rho=.0023769*exp ((-1*32.17*altft) / (1716* (oat+460) )); 
Seer—.5*rho* (vtrue%*2) ; 
const1=—-1* (qbar*sref) / (gw*vtrue/32.17); 
const2=qbar*sref*cbar/iyy; 
const3=const2*cbar/ (2*vtrue) ; 
const5=qbar*sref/gw; 


cbar=input (’AIRCRAFT MEAN CHORD LENGTH et 
CLa=input (‘INPUT DERIVATIVE CL a 1/RAD mee 
CMa=input (’ INPUT DERIVATIVE CM a 1/RAD or; 
CMg=input (’ INPUT DERIVATIVE CM q 1/RAD 4 
f ° 
f ° 


) 
) 
) ° 
) 
) 


™e Me Me Me 


lw=1750; % SCALE OF TURBULENCE 
menaltft<1750 
lw=altft; 
end 
sigw=08; % RMS VALUE OF TURBULENCE IN FT/SEC 
lamda= (vtrue/lw) ;k=(3*sigw%2) *vtrue/ (pi*lw) ;beta=vtrue/ (sqrt (3) *lw) 
% 
a=[constl*CLa 1 0 const1*CLa O; 
SemstZz2*CMa const3*CMg 0 constZ2*CMa oF 
0 1 0 0 OQ; 
0 0 0 0 ile 
0 0 0 -(lamda*2) —-2*lamda]j; 
% 
b={[const1*CLde 6) G:; 
const2*CMde 0 ON 
0 0 0; 
0 sqrt (k) /vtrue OF 
0 (sqrt (k) /vtrue) * (beta-2*lamda) O]; 
% 
c=[1 0 0 0 0; 
0 1 0 0 OG; 
0 0 i 0 0; 
const5*CLa 0 0 0 O); 
% 
da=[0 0 O; 
0 0 0; 
0 0 0; 
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const5*CLde 0 alge 


% 

&%SIMDATA(:,2)=AOA (RAD) SIMDATA (:,3)=PITCH RATE (RAD/SEC) 

%SIMDATA(:,4)=THETA (RAD) SIMDATA(:,5)=NORMAL ACC (G) 

% RAND NUMBER GENERATOR MEAN=0 VARIANCE=1 
rand(’normal’); 

{[phi, gam] =c2d(a,b,dt) ;u=[simdata(:,1),sysn*rand(ndp,1),ones(ndp,1) J 

[ynoise, xnoise]=dlsim(phi,gam,c,d,u) ; 

% STATENOISE + MEASUREMENT NOISE 

rand (” unitorm ))r- 

simdata(:,2:3:4:5)=ynoi1se:, 1: 2:2340e. 

+sysn*[.005618* (rand(ndp,1)-— 3s) 02903" ranegep ae) 

.005818* (rand(ndp, 1)=25) .0JS(randindpa i see 

% PLOTS FOR VIEWING ON MONITOR AND STORED IN META FILE 

*% ELVATOR vs TIME 

plot (t, (180/pi) *simdata(:,1)); 

xlabel(’Time (seconds)’);ylabel(’Elevator Input (degrees) ’) ; 

ans=[‘'title(’’’,typac,’ SIMULATED INPUT’’);']J; 

eval (ans) ;pause 

meta A:\plots\deltae 

% AOA vs TIME 

plot (t, (180/pi) *simdata(:,2)); 

xlabel(’Time (seconds)’);ylabel(’AOA Output (degrees) ’); 

ans=[‘’title(’’’,typac,’ SIMULATED DATA’’);']; 

eval (ans) ;pause 

meta A:\plots\AOA 

% PITCH RATE (Q) vs TIME 

plot (t, (180/p1) *simdatatcs))- 

xlabel(’Time (seconds)’);ylabel(’Pitch Rate, Q, Output (deg/sec)’); 


ans=[‘’title(’’’,typac,’ SIMULATED DATA’’);']J; 

eval (ans) ;pause; 

meta A:\plots\Q 

% THETA vs TIME 

plot (t, (180/pi) *simdata(:,4)); 

xlabel(’Time (seconds)’);ylabel(’Theta Output (deg)’); 
ans=['’title(’’’,typac,’ SIMULATED DATA’’);']J; 

eval (ans) ;pause; 

meta A:\plots\theta 

% NORMAL ACC vs TIME 

plot (t, SimdatatG, 5) > 

xlabel(’Time (seconds)’) ;ylabel(’Normal Acceration Output (G)’); 
ans=[’title(’’’,typac,’ SIMULATED DATA’’);']jJ; 

eval (ans) ;pause; 

meta A:\plots\G 

disp(’ % 

isp (0 p33 9 8 ee ee 
disp ( we 

disp (’NOTE ====> DATA BEING SAVED TO A .mat FILB’ ) 
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disp (' FOR MMLE PROCESSING ') 
if sysn>0 
ans=[‘’save N’,typac,’ typac simdata sysn truth iyy gw sref cbar’]; 
eval (ans) 
disp(’File name is N followed by the type aircraft.mat N’),typac; 
else 
ans=[’save NN’,typac,’ typac simdata sysn truth liyy gw sref 
char’]; 
eval (ans) 
disp(’File name is NN followed by the type aircraft.mat 
NN), typac; 
end 
Gusp(’ ‘) 
'dir/w 
amaat’ ! 
disp ° 
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Be LATDIR.M 


clear; 
*% MACRO FILE NAME >======= LATDIR.M =======< 


disp (ans) 

disp (74) 

disp(’ MACRO TO GENERATE SIMULATED LATERAL-DIRECTIONAL DATA’ ) 
disp (’ USING AILERON & RUDDER INPUTS WITH OR W/O NOISE’ ) 
disp (ans) 

disp(’ ‘) 

§ = — — - — — GET AIRCRAFT TYPE TO BE USED ------------ 

disp (’ AIRCRAFT TYPES AVAILABLE ‘s) 
dwapi(*) 3) 

disp (’NAVION A4D FI104A JETSTAR B747 UAV OTHER’ ) 
disp(’ ’) 

disp (‘NOTE ========> PROGRAM IS CASE SENSITIVE <========' ) 
disp(’ ’) 

disp(’SELECT "OTHER" TO INPUT DATA FOR USER DEFINED AIRCRAFT’ ) 
qisp( fa 


typac=input (’TYPE IN DESIRED A/C FROM THE ABOVE LIST. ', sar 
$ poo DETERMINE IF NOISE IS WANTED 

disp(’ ‘) 

sysn=input (’INPUT A 1 TO INCLUDE NOISE AND A ZERO FOR NO NOISE. ’); 
ndp=301;%------ NUMBER OF DATA POINTS - 15 SEC 

dat=.05; $------= TIME STEP FOR THE DATA 

amp=.05; $--—-=— AMPLITUDE OF INPUT (RADIANS) 

period=1; %----- PERIOD OF DOUBLET IN SEC = PERIOD + 1 
t=[0:ndp-1] *dt; $------------- TIME VECTOR 
Simdata=zeros (ndp, 7) ;%------- SETUP DATA MATRIX ALL ZEROS 

$n rrr GENERATE THE INPUT DOUBLET 


dslope=(4*amp*dt) /period; 

dl=(0:-1*dslope:-1l*amp) ;dla=-amp*ones (1:10) ; 
d2=(-1*amp+dslope:dslope: 0); 

d3= (dslope:dslope:amp) ;d3a=—-1*dla; 
d4=(amp—-dslope:-1*dslope: 0) ; 

§ wenn nnn nnn enn AILERON AND RUDDER INPUTS 
Simdata(:,1)=[zeros(1,60) dl dla d2 d3 d3a d4 

zeros (1,ndp-—(period/dt) —-81)]’; 

Simdata(:,2)=[-dl -dla -d2 -d3 -d3a -d4 

zeros (1, ndp—(pemrod, dty—21)) 


vtrue=input (’ INPUT AIRCRAFT TRUE VELOCITY ft/sec ‘Ve 
altft=input (’ INPUT AIRCRAFT ALTITUDE fe ‘es 
oat=input (‘INPUT THE OUTSIDE AIR TEMPERATURE Deg F "); 
if (stremp (’ CTHER’ 7 typace v7, 
typac=input (’ INPUT THE A/C TYPE ....< 6 characters ', soe 
sref=input (’ INPUT THE AIRCRAFT REFERENCE AREA ft%*2 
gw=input (‘INPUT THE AIRCRAFT GROSS WEIGHT lbs '); 


ixx=input (‘AIRCRAFT Ixx MOMENT OF INERTIA slug-£t*2 “se 
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ixz=input (’AIRCRAFT Ixz MOMENT OF INERTIA 
1Z2Z=input (’AIRCRAFT Izz MOMENT OF INERTIA 
bbar=input (’AIRCRAFT WING SPAN LENGTH "b" 


Saugene 2 
Sivug sft «2 


~ ~ 


~ 


% 
CYb=input (‘INPUT DERIVATIVE CY_b 1/RAD ae 
Clb=input (‘ INPUT DERIVATIVE Cl b 1/RAD ')>; 
CNb=input (‘’ INPUT DERIVATIVE CN b 1/RAD '); 
Clp=input (’ INPUT DERIVATIVE Cl _p 1/RAD ae 
CNp=input (’ INPUT DERIVATIVE CN p 1/RAD "); 
Clr=input (‘INPUT DERIVATIVE Cl r 1/RAD '); 
CNr=input (’ INPUT DERIVATIVE CN r 1/RAD ry); 
Clda=input (’ INPUT CONTROL DERIVATIVE Cl da 1/RAD / ee 
CNda=input (’ INPUT CONTROL DERIVATIVE CN da 1/RAD aya 
CYdr=input (‘INPUT CONTROL DERIVATIVE CY dr 1/RAD ae 
Cldr=input (‘INPUT CONTROL DERIVATIVE Cl dr 1/RAD oye 
CNdr=input (’ INPUT CONTROL DERIVATIVE CN dr 1/RAD "); 
% 
ans=[’save ’,typac,’ .mat;’ J; 
eval (ans) ;%---------- SAVE THE NEW A/C DATA IN A .MAT FILE 


mamieti— | cr Clb CNbaClp CNp Clr CNr Clda CNda CYdr Cldr CNdr]; 
else 

ans={’ load '’,typac,’ .mat;’ J; 

aa atis ) ;>——————— —— LOADS THE ‘TRUTH DATA’ ON SELECTD A/C 

meaen—|~CiYb Clb CNb Clp CNp Clr CNr Clda CNda CYdr Cldr CNdr]; 
end 
ae ee ee ee ee me ee ee ce ee CALCULATE DENSITY, DYNAMIC PRESSURE, AND MASS 
emo—.0023769*exp( (-1*32.17*altft) /(1716* (oat+460))); 
gqbar=.5*rho*vtrue*vtrue; 
mass=gw/32.17; 
a — CALCULATE CONSTANTS 
const1=(qbar*sref) /mass; 
const2=qgbar*sref*bbar;const2a=const2*bbar/ (2*vtrue) ; 
const3=const1/32.17; 


Sa el eee DRYDEN TURBULENT MODEL CONSTANTS 
mea ltftt<1750 
lw=altft; 
else 
tw=1750 ; $-------- SCALE OF TURBULENCE 
end 
aaa— 02; 5——————————— RMS VALUE OF TURBULENCE IN FT/SEC 


lamda=(vtrue/1lw) ;k=(3*sigw%2) *vtrue/ (pi*lw) ;beta=vtrue/ (sqrt (3) *lw) 


% SYSTEM STATE SPACE MATRICES 
% INERTIAL MATRIX 
if 


n=[1 0 CO ee Ort: 0% 
Oerxesx —1xz 0 O OQ; 
Oe=1x2 122 «0 0 =O; 
0 0 Gm is CO: . 0; 
0 0 Cy Ors 0; 
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0 0 Om 0 OO Ae 
an=[(const1*Cyb/vtrue) 0 -1 £=(32.17/vtrue) (const1*CYb/vtrue) 
0; 

(const2*Cl1lb) (constZ2a*Clp) (const2a*Clr) 0 (const2*Clb) 


0; 
(const2*CNb) (constZa*CNp) (constZa*CNr) Q £(const2*CNb) 
Os 
0 1 0 0 0 
0; 
0 0 0 0 0 1 
0 0 0 0 - (lamda%*2) 
-2*lamda]; 
% 
bn=(0 (const1*CYdr/vtrue) 0 Or 
(const2*Clida) (const2*Cldr) 0 OF 
(const2*CNda) (const2*CNdr) 0 OF 
0 0 0 UF 
0 0 (sqrt (k) /vtrue) O; 
0 O (sqrt (k) /vtrue) * (beta-2*lamda) Ol]; 
% 
a=inv (In) *an;b=inv (In) *bn; 
c=([(1 0 0 0 0 0; 
0 al 0 0 0 OF 
0 0 alt 0 0 8) 
0 0 0 il 0 0; 
const s*€Ybp ee 0 0 0 Gals; 
% 
d=[0 0 0 Oe 
0 0 0 Os 
0 0 0 0; 
0 0 0 0: 
0 const3*CYdr 0 Olt 
% 
% SIMDATA(:,1)=AILERON INPUT (RAD) SIMDATA(:,2)=RUDDER INPUT (RAD) 
% SIMDATA(:,3)=BETA (RAD) SIMDATA(:,4)=ROLL RATE (RAD/SEC) 
% SIMDATA(:,5)=YAW RATE (RAD/SEC) SIMDATA(:,6)=ROLL ANGLE (RAD) 
% SIMDATA(:, 7)=LATERAL ACC 
(phi, gam) —C2da, > ee) 0a CONVERT TO DISCRETE 
§-<—<--------- RAND NUMBER GENERATOR MEAN=0 VARIANCE=1 


rand(’normal’ ); 
u=(simdata(:,1),simdata(:,2),sysn*rand(ndp,1),ones(ndp,1)]; 
%--—* INPUTS 
(OUTY OUTX]=dlsim(phi,gam,c,d,u); 
Simdata (773-7) =OUTRY(. 4 5); 
rand(’uniform’ ); 
eee ADD THE MEASUREMENT NOISE 
Simdata(:,3:7)=simdata(:,3:7)+ 
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arn Oe soe rand (ndp,1)—.5))202909* (rand (ndp,1)-.5) 

Sees ~ (rand(ndp,1)—.5) .005818* (rand (ndp,1)-.5) 

fee rand (ndp,1)=—.5) ]; 
% PLOTS FOR VIEWING ON MONITOR AND STORED IN META FILE 
% AILERON INPUT 
Smeplot (211) plot (t, (180/pi) *simdata(:,1)); 
xlabel(’Time (seconds)’);ylabel(’Aileron Input (degrees)’); 
ans=[’title(’’’,typac,’ SIMULATED INPUT’’);’);eval(ans); 
% RUDDER INPUT 
Slee lot (212) ;plot (t, (180/pi) *simdata(:,2)); 
xlabel(’Time (seconds)’);ylabel(’Rudder Input (degrees) ’); 
ans=[’title(’’’,typac,’ SIMULATED DATA’’);’];eval (ans); 
pause; %meta A:\plots\Latinput 
% SIDE SLIP (Beta) 
eeprot (111) ;plot (t, (180/pi) *simdata(:,3)); 
xlabel(’Time (seconds)’);ylabel(’Side Slip ,B, Output (deg) ’); 
ans=[/’title(’’’,typac,’ SIMULATED DATA’’);’];eval (ans); 
pause; %meta A:\plots\sslipout 
% ROLL RATE (P) 
meet, (180/pi) *simdata(:,4)); 
xlabel(’Time (seconds)’);ylabel(’Roll Rate ,P, Output (deg/sec)’); 
ans=[’title(’’’,typac,’ SIMULATED DATA’’);’J];eval (ans) ; 
pause; %*meta A:\plots\rollrout 
% YAW RATE (R) 
Pleert, (180/pi) *simdata(:,5)); 
xlabel(’Time (seconds)’);ylabel(’Yaw Rate ,R, Output (deg/sec)’); 
mee citle(’’’,typac,’ SIMULATED DATA’’);’];eval (ans); 
pause; %meta A:\plots\yawout 
% ROLL ANGLE (phi) 
plot (t, (180/pi) *simdata(:,6)); 
xlabel(’Time (seconds)’);ylabel(’Roll Angle Output (deg)’); 
ans=[’title(’’’,typac,’ SIMULATED DATA’’);’]j;eval (ans) ; 
pause; meta A:\plots\rollaout 
% LATERAL ACCERATION (G) 
plot (t,simdata(:,7)); 
xlabel(’Time (seconds)’);ylabel(’Lateral Acceleration (G)’); 
ans=[/’title(’’’,typac,’ SIMULATED DATA’’);’)]);eval (ans) ; 
pause; %meta A:\plots\gout 
elc 
ete a nee Cu, 
esp (’ ’) 
disp (’NOTE ====> DATA BEING SAVED TO A .mat FILE’ ) 
ars (’ FOR MMLE PROCESSING ’) 

if sysn>0 

ans=[’save N’,typac,’ typac simdata sysn truth ixx 12z 1xz gw sref 
bbar’ ]; 

eval (ans) 
else 

ans=[’save NN’,typac,’ typac simdata sysn truth ixx 12z ixz gw 
sref bbar’]; 

eval (ans) 
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end 

disp a=) 

!'dir/w 

disp(’ NOW RUN npsmmle2.m WITH THE CREATED DATA FILE.’ ) 


68 


Cc. SIMULATED LONGITUDINAL MMLE 


1. NPSMMLE1.M 


% MACRO NAME >====== NPSMMLE!1.M =======< 
% Date: 3 Feb 92 
clear; 


ferase npsmmlel.log; 

diary npsmmlel.log 

a 1s 

disp(’ NPS PARAMETER IDENTIFICATION MACRO FILE ’) 

disp(’FOR SIMULATED FLIGHT DATA USING SIMPLFIED SHORT PERIOD ’) 


disp (’ LONGITUDINAL STABILITY AND CONTROL DERIVATIVES’ ) 
ausp(’ ’) 

ea ———————————-——- — ——- —- —- — = = nn es 
mogmmat ) $—-———————==--------— RUN INITIALIZATION MACRO 


format compact,clec 
load npsinitl; 
global sref cbar gw iyy vtrue qbar dt all ql thi anl del; 


nn rrr rr rrr rrrrc--- INPUT FLIGHT TEST DATA 
Wveata—zeros (ndp, 6) ;% -------- ESTABLISH DATA MATRIX FOR MMLE.M 
$a nn nnn nnn nnn UYDATA(:,1) = DELTA E (INPUT) 
Gn rrr rrr rrrsrsss UYDATA(:,3) = AOA 

== —— ——— —— — ———— —- —-—--_---------- UYDATA(:,4) = PITCH RATE (Q) 
a a ea UYDATA(:,5) = THETA 

ee eee UYDATA(:,6) = NORMAL ACC 
=< —— COLUMN NUMBER ONE IN DATA FILE ELEVATOR INPUT 
coll=[’uydata(:,1)=simdata(:,1);’]J;eval(coll); 

Eee COLUMN NUMBER TWO IN DATA FILE UNITY INPUT 
uydata(:,2)=ones (ndp,1); 

——--—— COLUMN NUMBER THREE IN DATA FILE ANGLE OF ATTACK 
col3=[’uydata(:,3)=simdata(:,2);'’];eval(col3) ; 

§e——--— COLUMN NUMBER FOUR IN DATA FILE PITCH RATE 
col4=[’uydata(:,4)=simdata(:,3);'’];eval(col4) ; 

§---- = COLUMN NUMBER FIVE IN DATA FILE THETA 
col5=[’uydata(:,5)=simdata(:,4);’];eval (col5) ; 

§ee---- COLUMN NUMBER SIX IN DATA FILE NORMAL ACC 
col6=[’uydata(:,6)=simdata(:,5);’];eval(col6) ; 

§<------- INITIAL CONDITIONS 


del=uydata(1,1);all=uydata (1,3) ;ql=uydata(1,4); 
thl=uydata (1,5) ;anl=uydata (1, 6); 
——————— ADDITIONAL INPUTS TO MMLE FOLLOW 
mesmam='npsp2ssl’; * -------- MACRO NAME FOR P2SS FUNCTION 
pO=pref; %*- INITIAL PARAMETER ESTIMATES INPUT DURING NPSINIT1.M 
%$-- CHECK IF THE WEIGHTING FUNCTION IS TO BE USED FOR INITIAL 
VALUES 
disp(’DO YOU WANT TO WEIGHT THE INITIAL ESTIMATES? INPUT 1=Y O=NO 
f 

) 
input (' '); 
if ans== 
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disp(’Input Weighting Row Vector length 1 x 5 ’) 
disp(’Use brackets- ex. [.1 11 .11] & lower # higher weight’ ) 
rmsOQ=input (’ ’); 


end 
pidq=[1];%--- IDENTIFY WHICH PARAMETERS ARE TO BE IDENTIFIED 
pidm=[1:5];%--- IN THE QUADRATIC, MARQUARDT, AND FINAL STAGES. 


pidf=[1:5];%- pidq, pidm, pidf MUST BE VALID EVEN IF NOT USED 
opt=([0 5 10 10 (02 (00S 00 ei, 

if sysn==0, opt (4) =0;end 

$- DEFAULT ITERATIONS AND CONVERGENCE CRITERIA, IF NOISE FREE 
OPT4=0 

gg0=eye (4) *(.001) ; $------------- INNOVATIONS COVARIANCE MATRIX 
pert=le-—4;%-PERTURBATION USED FOR NUMERICAL GRADIENT CALCULATION 
linesearch=1;%--- USE LINESEARCH TO HELP PROC. NOISE CONVERGENCE 
!cls 

mmle $-o-oo torr CALL MAIN MMLE MACRO FROM TOOL BOX 
ee PERFORM FINAL CALCULATIONS 
cla=pfin(1);cma=pfin (2) ; cmq=pfin (3) ; clde=pfin (4) ; cmde=pfin (5); 
disp(’ '’) 

deriv=[cla cma cmq clde cmde]; 


disp (’ MMLE STABILITY & CONTROL DERIVATIVES | 
disp (’ CLA CMA CMQ CLDE CMDE’” ) 
disp (deriv) 

disp (’ INITIAL INPUT DERIVATIVES =) 


disp (pref) 

if exist (’truth’ ) ==1; 
disp(’ TRUTH DERIVATIVES USED TO GENERATE THE DATA ’); 
disp (truth) 

end 

pause 

mleplotl 

diary off 

$! print npsmmlel.log; 

+{————————___ _ _ _ _ — __ ___ soci a END NPSMMLE1 .M 
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2. NPSINIT1.M 


%& MACRO FILE NAME >=2===== NPSINIT1.M =======< 

% Date: 3 Feb 92 

% INITIAL SETUP MACRO FOR RUNNING THE PARAMETER 

% ESTIMATION TOOLBOX IN MATLAB FOR SIMPLIFIED LONGITUDINAL 

% SHORT PERIOD STABILITY AND CONTROL DERIVATIVES 

% THIS MACRO ‘GETS’ CONSTANTS AND DATA FILE 

§-<-——— LOAD DATA FILE 

disp(’ DATA FILE MUST CONTAIN DATA IN COLUMN MATRIX '’) 

disp(’ MATRIX NAME SIMDATA: N DATA PTS X 5 COLUMNS ’) 

disp(’ ELEVATOR/ALPHA/THETA/PITCH RATE/NORMAL ACC.’) 

disp(’ “’) 

disp(’ DATA FILE NAME--MUST EXIST WITH A .mat EXTENSION. ’); 
data=input (’ENTER DATA FILE NAME ( ==> WITH OUT <== .MAT 
EXTENSION)? /’,'s"); 

if exist (’dt’)==0,dt=input (’DELTA T BETWEEN DATA POINTS? ’);end 
ldc=[’load ’,data,’.mat;’]; 

eval (ldc) ;$--EXECUTES LOAD COMMAND 

ndp=length (simdata) ; 

See ndp—1] *dt ; 

§ eee INPUT REQUIRED CONSTANTS--—- IF NOT IN DATA FILE 

if exist (’sref’)==0,sref=input (’REFERENCE AREA (S) IN SQUARE FEET? 
‘yY;end 

if exist (’cbar’ )==0, cbar=input (’MEAN AERODYNAMIC CHORD IN FEET? 
‘yY;send 

if exist ('gw’ )==0,gw=input (‘AIRCRAFT GROSS WEIGHT IN POUNDS? ’);end 


if exist (/iyy’)==0,iyy=input (‘MOMENT OF INERTIA (IYY) IN SLUG-FT*2? 
‘y;end 
$sdc=[’save ’',data,’ simdata dt sref cbar gw iyy ’'];eval(sdc); 


vtrue=input (’AIRSPEED IN FEET PER SECOND? '); 
altft=input (‘AIRCRAFT ALTITUDE IN FEET? ‘'); 

oat=input (’OUTSIDE AIR TEMPERATURE? Cie 

$------ CALCULATE CONSTANTS DENSITY AND DYNAMIC PRESSURE 
rho=.0023769*exp ( (-1*32.174*altft) / (1716* (oat +460) )); 
qbar=.5*rho*vtrue*vtrue; 

$--------------- INPUT INITIAL ESTIMATES FOR STABILITY 
%$--------------- AND CONTROL DERIVATIVES 

pref (1)=input (‘CL ALPHA ESTIMATE (1/RAD)? '); 

pref (2)=input (‘CM ALPHA ESTIMATE (1/RAD)? ’); 

pref (3)=input (‘CM Q ESTIMATE (1/RAD)? ’) 

pref (4)=input (’CL DE ESTIMATE (1/RAD)? ’ 
pref (5)=input (‘CM DE ESTIMATE (1/RAD)? ’ 
l!erase npsSinitl.mat; 

Save npsinitl 
%------------------------------------------ END NPSINIT1.M 


f 


f 


) 
) 


yal 


3. NPSP2SS1.M 


function [a,phi, gampe7a,a,.t.ae pn oe ig at 


% MACRO FILE NAME >==S=S>=== NPSP2SS1. M ======>=< 

% Date: 3 Feb 92 

%; SS Ee 
% MACRO TO EST. FUNCTION FOR TRANSFORMING 

% MODEL PARAMETERS INTO STATE SPACE EQUATIONS 

% NN SS 
%§ P2SS FUNCTION FOR NPSMMLE1.M 

5 ee 
% p(l) = CL ALPHA STABILITY AND CONTROL 

$ p(2) = CM ALPHA PARAMETERS 

% p(3) = CM Q | 

% p(4) = CL DE | 

SP (5)  =eGMEbE | 

$m nnn nn nnn eee PERFORM INITIAL CALCULATIONS 


const1=(-1*qbar*sref) / (cos (all) *gw*vtrue/32.17) ; 
const2=qbar*sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 

const 4=32.17%*cos (thl1) / (vtrue*cos (all) ); 
const5=(qbar*sref) /gw; 

const 6=-1*% (const1*p (1) *all+qi+constl1*p (4) *del+const4) ; 
const 7=-1* (const2*p (2) *all+const3*p (3) *qitconst2*p (5) *del) ; 
const@= (-1* (const5*p (1) *all+const5*p (4) *dem™) +ani; 


$----------------------- -- --- -------- STABILITY DERIVATIVES 
a=([const1*p (1) il 0; 

COnSEZ*p (Zz) const 3*pi3) O; 

0 1 Oj; 
$------------------------- ----------- CONTROL DERIVATIVES 
b=[(const1*p (4) (const4+const6) ; 

const2*p (5) const7; 

Q (=1*aaeie, 

a a aa MEASUREMENT MATRIX 
e—1 1) 0 0; 

0 ie OF 

0 0 a 

COnStES~pi 1) 0 Oj; 
fern nn nr rn nr nn nn nnn rr rer sr serco=-— FEED THROUGH MATRIX 
dad=[0 0; 

0 0; 

0 O; 

const5*p (4) const8]; 
%-------------------------------------- STATE NOISE COVARIANCE 
g=eye (a) *le-4;%-- Q IS THE SAME SIZE AS a 
%---------------- WITH Q*Q’ POS. DEFINITE 
%---------------- ROWS IN Q IN WHICH PARAMETERS OCCUR, A VECTOR 
% SAME DIMENSION AS p 


rowing=[(0 0 0 O Oj; 
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aaa. INITIAL STATE VECTOR 


%—-—--—------------- DISCRETIZE 

% *x****NEED TO EDIT dt (below) FOR THE DELTA T OF THE DATA ***** 
a@e—.05; 

[phi, gam]=c2d(a,b,dt); 

0 See SS Se Se SS SS a END NPSP2SS1.M 


a3 


4. MLEPLOT1.M 


% MACRO FILE NAME >======= MLEPLOTI1.M =======< 
% Date: 3 Feb 92 
% -—---- MACRO TO PLOT DATA FROM NPSMMLE1 


const1l=(-1*qbar*sref) / (cos (all) *gw*vtrue/32.17) ; 
const2=qbar*sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 
const4=32.17*cos (th1) / (vtrue*cos(all)); 

const 5=(qbar*sref) /gw; 

const 6=-1* (const1*pfin (1) *all+ql+const1*pfin (4) *del+const4) ; 

const 7=-1* (const2*pfin (2) *all+const3*pfin (3) *qi+const2*ptintS aoa 


const 8= (-1* (const 5*pfin (1) *all+const5*pfin (4) *del))+anl1; 


:-Soetesteetasteeateesteeatestanetenetoeteesteeatane tenet eetoestanleectease teeta teertetenten entation STABILITY DERIVATIVES 
a=[const1*pfin (1) 1 OF 
GOnstZ-prrn) Const34 pian (>) O; 
0 al Ol; 
5 ee a a = = SS = CONTROL DERIVATIVES 
b=[const1*pfin (4) (const4+const6) ; 
const2*pfin (5) const7; 
0 (-1*ql) J; 
$m nnn MEASUREMENT MATRIX 
c=[1 0 CG; 
0 1 0; 
0 0 Ls 
CONSt>*<pLrnit |) 0 Oj): 
Seana ehh anette ieee etait anette ete teententeneenetete FEED THROUGH MATRIX 
d=[0 OF 
0 OQ; 
0 Or 
Const 5*peEin (4) const8]j; 
rtdc=(ils07 pil); 
$e eee ee OUT2 = OUTPUT VECTOR---- OUT3 = STATE VECTOR 


[OUTZ, OUT3)=lsim(a,b,c,d,uydata(:, 1:72), e407), 

if exist (’truth’ )>0; 

const 6=-1* (const1l*truth (1) *all+ql+constl*truth (4) *del+const4) ; 
const 7=-1* (const2*truth (2) *all+const3*truth (3) *ql+constZ*trumie a 
el); 

const 8= (-1* (const 5*truth (1) x*all+const5*truth (4) *del))+anl1; 


a=[constl*truth (1) 1. OF 
const2*truth (2) const3*truth (3) OF 
0 1 O]; 
SS aaa la eee CONTROL DERIVATIVES 
b=[constl1*truth (4) (const4+const6) ; 
const2*truth (5) const7; 
0 (legal) || ¢ 
$—--—--—--— nnn MEASUREMENT MATRIX 
c=[1 0 O- 
0 i QO; 
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0 0 ILS 


const5*truth (1) 0 Oo); 
%------------------------------------- FEED THROUGH MATRIX 
d=[0 0; 

0 0; 

0 0; 

const5*truth (4) const8]; 

fru2, TRU3)=lsim(a,b,c,d,uydata(:,1:2),t,x0); 
end 
5 PLOTS TO MONITOR AND STORED IN META FILE 


if exist (’typac’)==0,typac=input (’ INPUT THE AIRCRAFT TYPE ? 
wes ) send 
ferase a:\plots\outputg.met; 


lerase a:\plots\outputth.met; 
!'erase a:\plots\outputde.met; 
f'erase a:\plots\outaoa.met; 


ferase a:\plots\outputq.met; 
$----ELEVATOR VS TIME 
Wemewort, plot (t, rtdc*uydata(:,1),'~—r’); 
xlabel(’Time (seconds)’);ylabel(’Elevator Input (degrees)’); 
ans=[/title(’’’,typac,’ ELEVATOR INPUT VS TIME '’);']; 
eval (ans) ;pause 
meta A:\plots\OUTPUTde 
§—-—-——— AOA (OBSERVED AND ESTIMATED) VS TIME 
Parorsee, PCAc*uydata(:,3),’*r’');hold on; 
xlabel(’Time (seconds)’);ylabel(’AOA (degrees) ’) ; 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED AOA RESPONSE’’);']; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’,’sc’);pause; 
Pome, TCAC*OUTZ (:,1),’ og’ ); 
text (.6,.80,’o Estimated Response ’,’sc’) ;pause; 
gmeesxist (’ truth’ )>0; 
Pact, TLAOC*TRUZ(:,1),'-b’); 
text (.6,.75,’- True Response ’,’sc’) ;pause; 
end 
pause 
meta A:\plots\outAOA 
5-—--—-— QO (OBSERVED AND ESTIMATED) VS TIME 
Wemomort -plot (t, rtdce*uydata(:,4),'*r’);hold on; 
xlabel(’Time (seconds)’);ylabel(’Pitch Rate, Q, (deg/sec)’); 
ans=[/title(’’’,typac,’ ESTIMATED AND MEASURED q RESPONSE’’);']; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’,’S8sc’);pause; 
peo, rCAc*OUTZ (:,2),' og’); 
text (.6,.80,’o0 Estimated Response ’,’sc’);pause; 
memexist (‘’truth’ )>0; 
Meee (tc, LEQC*TRUZ(:,2);,'’-—b’); 
text (.6,.75,’- True Response ’,’sc’);pause; 
end 
pause 
meta A:\plots\OUTPUTQ 


Jas, 


4—--—=-= THEATA (OBSERVED AND ESTIMATED) VS TIME 
hold off; plot (t, rtde*uydata(:,5))) “© )7 moter eu, 
xlabel(’Time (seconds)’);ylabel(’Pitch Angle, Theta, (deg)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED THETA 
RESPONSE’ ’);' J; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’,’sc’);pause; 
plot (t, rtde*OUT2Z (2735) coe. 
text(.6,.80,’o Estimated Response ’,’sc’) ;pause; 
if exist (’ truth eo 

plot (€, rtde* TRU ss), a, 

text(.6,.75,’- True Response ’,’sc’) ;pause; 
end 
pause 
meta A:\plots\OUTPUTTH 
j= ACCELERATION (OBSERVED AND ESTIMATED) VS TIME 
hold off;plot(t, uydata(:, 6), “r" )) notaen, 
xlabel(’Time (seconds)’);ylabel(’Acceleration, G ’); 
ans=['title(’’’,typac,’ ESTIMATED AND MEASURED G RESPONSE’’);']; 
eval (ans) 
text(.6,.85,’'* Measured Data Points ’,’sc’) ;pause; 
PploticC, OUTZ (ay 4) cece). 
text (.6,.80,’o Estimated Response ’,'’sc’) ;pause; 
tf exist ( Eruen ) 0; 

plot (©, BRUCK: 4 bee 


text (.6,.75,’'- True Response ’,’sc’);pause; 
end 
pause 
meta A:\plots\OUTPUTG 
held off - 
Sada ee hh ae aa eee ea END MLEPLOT1.M 
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D. SIMULATED LATERAL-DIRECTIONAL MMLE 


1. NPSMMLE2.M 


% MACRO NAME >====== NPSMMLEZ .M =222=22==< 
% Date: 3 Feb 92 
clear; 


ferase npsmmle2.log; 


disp (ans) 

disp(’NPS STABILITY PARAMETER IDENTIFICATION MACRO FILE’ ) 
disp (’ FOR LATERAL-DIRECTIONAL 2) 
asp (' STABILITY AND CONTROL DERIVATIVES .) 
Gasp ('«‘) 

disp (ans) 

npsinit2 %----- RUN INITIALIZATION MACRO 


hermat Compact, clc; 

load npsinit2; 

global sref bbar gw ixx ixz izz vtrue qbar dt betal rolll yawl 
Bamgtel ayl dal dri; 


———— — — — — — - -- --- INPUT FLIGHT TEST DATA 
uydata=zeros (ndp, 8) ;% --- ESTABLISH DATA MATRIX FOR MMLE.M 

% UYDATA(:,1) = DELTA Aileron (INPUT) UYDATA(:,5) = ROLL RATE 
(p) 

% UYDATA(:,2) = DELTA Rudder (INPUT) UYDATA(:,6) = YAW RATE (r) 
% UYDATA(:,3) = UNITY INPUT UYDATA(:,7) = ROLL ANGLE 
(phi) 

% UYDATA(:,4) = BETA (SIDE SLIP) UYDATA(:,8) = LATERAL G 
(ay) 

$-—-—--—- COLUMN NUMBER ONE IN DATA FILE AILERON INPUT 
coll=[’uydata(:,1)=simdata(:,1);'’];eval (coll); 

---- COLUMN NUMBER TWO IN DATA FILE RUDDER INPUT 
col2=[’uydata(:,2)=simdata(:,2);'’J;eval(col2); 

uydata(:,3)=ones (ndp,1);% UNITY INPUT 

aa COLUMN NUMBER FOUR IN DATA FILE BETA (SIDE SLIP) 
col4=[’uydata(:,4)=simdata(:,3);’J;eval(col4); 

%---- COLUMN NUMBER FIVE IN DATA FILE ROLL RATE (p) 
col5=[’uydata(:,5)=simdata(:,4);'’J];eval(col5); 

%---— COLUMN NUMBER SIX IN DATA FILE YAW RATE (r) 
col6é=[’uydata(:,6)=simdata(:,5);’];eval(col6) ; 

$---- COLUMN NUMBER SEVEN IN DATA FILE ROLL ANGLE (phi) 
col7=[’uydata(:,7)=simdata(:,6);’J;eval(col7); 

%----— COLUMN NUMBER EIGHT IN DATA FILE LATERAL G (ay) 

col8=[(’ uydata(:,8)=simdata(:,7);'];eval(col8); 

eee INITIAL CONDITIONS FOR da, dr, BETA, p, r, phi, ay 


betal=uydata (1,4) ;rolll=uydata (1,5) ;yawl=uydata (1, 6); 
ranglel=uydata (1,7) ;ayl=uydata (1, 8); 
dal=uydata (1,1) ;drl=uydata (1,2); 


at 


+———[se ADDITIONAL INPUTS TO MMLE FOLLOW 


p2snam='npsp2ss2’ ;%----- MACRO NAME FOR P2SS FUNCTION 
p0O=pref;%---- INITIAL PARAMETER ESTIMATES 

$rmsQ=--- IF USED IT IS THE WEIGHTING FUNCTION FOR INITIAL VALUES 
pidg=— EL. +-—----—-—--+ IDENTIFY WHICH PARAMETERS ARE TO BE IDENTIFIED 
piGmeieed 2) s=——-— oe IN THE QUADRATIC, MARQUARDT, AND FINAL 
STAGES. 

pidt=(1+12);t=.-—se pidq, pidm, pidf MUST BE VALID EVEN IF NOT USED 


opt=[0 5 5 10 .02 .05 .001 1)];%- DEFAULT ITERATIONS AND CONVERGENCE 
% CRITERIA, IF NOISE FREE OPT4=0 
if sysn==0, opt (4) =0;end 


ggQ=evye (5) *(.01) ; 4-34. aaa eee INNOVATIONS COVARIANCE MATRIX 
pert=le~-4;%-- PERTURBATION USED FOR NUMERICAL GRADIENT CALCULATION 
linesearch=1; $----- USE LINESBARCH TO HELP PROC. NOISE CONVERGENCE 
Nn les--——=———_—_— ee CALL MAIN MMLE MACRO FROM TOOL BOX 

$9 === PERFORM FINAL CALCULATIONS 


CYb=pfin (1) ;Clb=pfin (2) ; CNb=pfin (3) ;Clp=pfin (4) ; 
CNp=pfin (5) ;Clr=pfin (6) ;CNr=pfin(7) ;Clda=pfin (8) ; 
CNda=pfin (9) ;CY¥dr=pfin (10) ;Cldr=pfin (11) ;CNdr=pfin(12); 


deriv=[CYb CED CNb Crp CNp Crass 

CNr Clda CNda CN. Cidr CNdr]; 
disp (’ MMLE STABILITY & CONTROL DERIVATIVES ry 
disp? CxS Clib CN b Clp CN p Cla 7 
disp(deriv(l1,:)) 
disp(’CN r Clda CNda Gr Gole Cilar CNdr’ ) 
disp (deriv(2,:)) 
disp (’ INITIAL INPUT DERIVATIVES ‘) 


disp (pref (1:6)); 
disp (pref (7:12)); 
if exist (‘’truth’ )==1; 
disp (’ "TRUTH DERIVATIVES" USED TO GENERATE DATA ‘ ) 
disp (truth(1:6)); 
disp(trutht 7.22 
end 
pause;mleplot2 
diary off 
$!print npsmmle2.1log; 
$e ee ee END NPSMMLE2 .M 
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2. NPSINIT2.M 


clear; 

* MACRO FILE NAME >======= NPSINIT2.M =======< 

% Date: 3 Feb 92 

%* INITIAL SETUP MACRO FOR RUNNING THE PARAMETER 

% ESTIMATION TOOLBOX IN MATLAB FOR LATERAL—DIRECTIONAL 

*% STABILITY AND CONTROL DERIVATIVES 

* THIS MACRO ’GETS’ CONSTANTS AND DATA FILE 

darsia(’ ') 

disp (’ DATA FILE NAME--MUST EXIST WITH A -mat EXTENSION. dap 3 
data=input (’ ENTER DATA FILE NAME ( ==> WITH OUT <== .MAT 


EXTENSION)? ‘','8s'); 
ide=(* load ',data,’ .mat;’];eval(ldc); 
ndp=length (simdata) ; 


beg exist (’ dt’) ==0,dt=input (’ INPUT dt BETWEEN DATA POINTS. 7) send 
Pe ndp—1]) *dt; 

5=—-=——-—-------- INPUT REQUIRED CONSTANTS 

if exist (’sref’)==0,sref=input (’REFERENCE AREA (S) IN SQUARE FEET? 
‘)y;end 

if exist (’bbar’ )==0,bbar=input (’WINGSPAN IN FEET? ’');end 

1f exist (’ gw’ )==0, gw=input (’AIRCRAFT GROSS WEIGHT IN POUNDS? ’);end 


if exist (’1xx’)==0,1ixx=input (’MOMENT OF INERTIA (Ixx) IN SLUG-FT%2? 
‘y;end 

if exist (‘ixz’)==0,ixz=input (‘MOMENT OF INERTIA (Ixz) IN SLUG-FT%*2? 
‘);end 

if exist (’12zz’ )==0,1zz=input (’MOMENT OF INERTIA (Izz) IN SLUG-FT%2? 
‘);pend 

$sdc=[(’save ’,data,’ simdata dt sref bbar gw ixx 1xz 12z 

‘];eval (sdc) ; 

vtrue=input (’ AIRSPEED IN FEET PER SECOND? ’); 

altft=input (’AIRCRAFT ALTITUDE IN FEET? '); 

oat=input (’OUTSIDE AIR TEMPERATURE? ) 

rho=.0023769*exp ( (-1*32.174*altft) / (1716* (oat+460) )); 
G@ear—.>*rho*vtrue*vtrue; 


We 


ee INPUT INITIAL ESTIMATES FOR STABILITY 


pref (10)=.175;%input (’CY dr FROM WIND TUNNEL (1/RAD) ? 
pref (11)=.02;%input (‘Cl dr FROM WIND TUNNEL (1/RAD)? ’ 
pref (12)=-.075;%input(’CN dr FROM WIND TUNNEL (1/RAD)? ’); 
!'erase npsinit2.mat; ” 

save npsinit2 

%----------------------- -- --- --- ------------ END NPSINIT2.M 


$< AND CONTROL DERIVATIVES TO START MMLE PROGRAM 
pref (1)=-.6;% input (’CY beta FROM WIND TUNNEL (1/RAD)? ’); 
pref (2) =-.15;% input(’Cl beta FROM WIND TUNNEL (1/RAD)? ’); 
pref (3)=.20;% input (’CN beta FROM WIND TUNNEL (1/RAD)? '); 
pref (4) =-.35;% input (’Cl_p FROM WIND TUNNEL (1/RAD)? '); 
pref (5)=-.05;% input (’CN p FROM WIND TUNNEL (1/RAD)? ’); 
pref (6)=.15;% input (‘Cl r FROM WIND TUNNEL (1/RAD)? '); 
pref (7)=-.2;% input (’CN r FROM WIND TUNNEL (1/RAD)? ’); 
pref (8)=.05;% input(’Cl da FROM WIND TUNNEL (1/RAD)? ’); 
pref (9)=-.001;% input(’CN da FROM WIND TUNNEL (1/RAD)? ’ 
) 
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3. NPSP2SS2.M 


maineeson [a,phi,gam,c,d,q,x0,dt, Sees pee ees ‘P) 


% MACRO FILE NAME >===S==== NPSP2SS2. M =======< 

% Date: 3 Feb 92 

% a ee ee ee ee es 
% MACRO TO ESTABLISH FUNCTION FOR TRANSFORMING 

% MODEL PARAMETERS INTO STATE SPACE EQUATIONS 

%& a a 
% P2SS FUNCTION FOR NPSMMLEZ2.M 

a 
ope) = CY beta (| ew) =s EN oz 
* p(Z) = Cl _ beta |STABILITY AND CONTROL | p(8) = Cl_da 
%$ p(3) = CN beta | | p(9) = CN da 
+ p(4) = Clp | PARAMETERS | p(l0) = CY dr 
%$ p(5) = CN | (Poel = C lar 
$p(6) = Clr | pep tt2) =. Cl dar 
% 


ee PERFORM INITIAL CALCULATIONS 
constl=(qbar*sref) / (gw/32.17); 
const2=qbar*sref*bbar;const2a=const2*bbar/ (2*vtrue) ; 
const3=(qbar*sref) /gw; 

$a wn nnn eee ee ee INERTIAL MATRIX 

In=[1l1 0 on 0; 

O aixx -1xz O; 

We=2xX2Z 122 OO; 


0 0 6 
an ee eee aeons PLANT 
an=[const1*p(1)/vtrue 0 -1 CS 2 ere): 
eenst2*p (2) constZa*p (4) constZ2a*p (6) O; 
BonstZ2*p (3) constZa*p (5) const2a*p (7) OF 
0 al 0 Oilee 
% 
a=inv (In) *an; 
% 
bn=[0 const1*p (10) /vtrue oe 
Sonstz*p (8) Gonsez~p(l1) (Os 
eGonst2*p (9) const2*p (12) OF 
0 0 0); 
% 
b=inv (In) *bn; 
% 
c=[1 0 0 O; 
0 1 0 O; 
0 0 i OG 
0 0 0 ee 
eonst3s*p (1) 0 0 Or: 
% 
d=[0 0 betal; 
0 0 EOL 1 
0 0 yawl; 
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0 0 ranglel; 


0 const3*p (10) ayl]; 
% 
§ == -=— + -— a a eee STATE NOISE COVARIANCE 
q=eye (a) *le-4 ; 4-=sa2=—=— Q IS THE SAME SIZE AS a 
% WITH Q*Q’ POS. DEFINITE! 
5 ———=. _—— =e ee ROWS IN Q IN WHICH PARAMETERS OCCUR, A VECTOR 
% SAME DIMENSION AS p 
rowing=0*p; 
5 eee INITIAL STATE VECTOR 
x0=[{betal rolll yawl ranglel]; 
a DISCRETIZE 
Gite 
[phi, gam)=c2d(a,b,dt) ; 
$e rrr rr nn nr rrr ssfcccc--— END NPSP2SS2.M 
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4. MLEPLOT2 .M 


% MACRO FILE NAME >S====== MLEPLOT2Z .M =======< 

% Date: 3 Feb 92 
———-<-. MACRO TO PLOT DATA FROM NPSMMLEZ2 
a GENERATE THE PREDICTED DATA 
aa. CALCULATE STABILITY AND CONTROL MATRIX 


mass=gw/32.17;rtdc=(180/pi) ;const1=(qbar*sref) /mass; 
const2=qbar*sref*bbar; const2a=const2*bbar/ (2*vtrue) ; 


const3=qbar*sref/gw; 
% 
In=[1 0 0 G 

O a1xx -1xz O; 

O -1xz izz O; 

0 0 0 1); 


an=[const1*pfin (1) /vtrue 0 -—1 B27 /verue: 
eenstZ*pfin (2) const2a*pfin(4) const2a*pfin (6) Oi; 
eenst2*pfin (3) const2a*pfin(5) const2a*pfin (7) Or. 
0 1 0 On 
a=inv (In) *an; 
% 
bn— [0 const1l*pfin(10) /vtrue OF 
PenstZ*pfin (8) GOMs22 primi...) 0 
eGonstZ2*pfin (9) GeonstZ*pfan (12) Oe 
0 0 Oj; 
b=inv (In) *bn; 
% 
ee (1 0 0 Oe 
0 Al 0 Oe 
0 0 ut Oe 
0 0 0 lee 
const3*pfin (1) 0 0 Oj; 
% 
d=([0 0 betal; 
0 0 Hol 1 
0 0 yawl; 
0 0 ranglel; 
0 COonses“pfain (10) ayl]; 
% OUTZ = OUTPUTS ---- OUTS = STATE VECTOR 
[OUT2, OUT3]=lsim(a,b,c,d,uydata(:,1:3),t,x0); 
if exist (’truth’ )>0 
an=[constl*truth (1) /vtrue 0 -1 Zee We Ale ¢ 
const2*truth (2) const2Za*truth (4) const2a*truth (6) 
0; 
eonstZ2*truth (3) constZa*truth(5) constZa*truth (7) 
0; 
0 1 0 Ol]; 


a=inv (In) *an; 
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% 


bn=)[0 constl*truth (10) /vtrue 0; 
const2*truth (8) const2*truth (11) O= 
const2*truth (9) const2*truth (12) 0; 
0 0 OF 

b=inv (In) *bn; 

% 

c=[1 0 0 OQ; 

0 1 0 0; 
0 0 i OF 
0 0 0 ils 
const3*truth (1) 0 0 Oj; 

% 

d=([0 0 O; 

0 0 O- 

0 0 OF 

0 0 0; 

0 Const 3*truth (lo) Oj; 
% 


ETRUZ = OUTPUT TRU3 = STATE VECTOR 

[TRUZ2, TRU3]=lsim(a,b,c,d,uydata(:,1:3),t); 
end 

% PLOTS FOR VIEWING ON MONITOR AND STORE TO META FILE 
f'erase a:\plots\*.met 
subplot (211) ; plot (t, rede*uvearat. |); 
xlabel(’Time (seconds)’);ylabel(’Aileron Input (degrees) ’); 
ans=[(’taitle(’’’,typac,’ AILERON INPUT VS TIME ’'’);'’);eval (ans) 
subplot (212) ;plot(t, rede*uvdaeat zen 
xlabel(’Time (seconds)’);ylabel(’ Rudder Input (degrees) ’); 
ans=[(‘/title(’’’,typac,’ RUDDER INPUT VS TIME ’’);'’j)];eval (ans) 
pause; %meta A:\plots\INPUTDAR 
% Beta vs Time 
subplot (111) ;plot (t, rtde*uydata(:,4)7 *rn") 7 noldger, 
xlabel(’Time (seconds)’);ylabel(’Beta (degrees)’); 
ans=['’title(’’’,typac,’ ESTIMATED AND MEASURED Beta 
RESPONSE!’ ’);')]; 
eval (ans) 
text(.6,.85,’* Measured Data Points ’,’sc’) ;pause; 
plot (C; sede*OUTZ(:) 1) eG an. 
text(.6,.80,’o Estimated Response ’,’sc’) ;pause; 

if exist (erurcn =)-0- 

plot (t; rede=TRUZ(. 1) soe 
text (.6,.75,'-  "Tnue Response au Gace.) 

end;pause; %meta A:\plots\outbeta 
% Roll rate vs Time 
hold off; ploett, rtdc*uydaeatey oe ie nota en, 
xlabel (’Time (seconds)’);ylabel(’Roll Rate, p, (deg/sec)’); 
ans=[/title(’’’,typac,’ ESTIMATED AND MEASURED p RESPONSE’’);']; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’,’sc’) ;pause; 
plot (t; rEde*OU Ta. 2), 109 a, 
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text (.6,.80,’o Estimated Response ’,’sc’);pause; 
mm exist (’truth’ )>0; 
Plotit, rtdc“tRuzt: ;, 2) ,"—b_):; 
Eon O, 75, —s)' True Response” ”,’ sc’); 
end; pause; tmeta A:\plots\OUTP 
% Yaw Rate vs Time 
femanort ;plot (t, rtdc*uydata(:,6),*’*r’) shold on; 
xlabel(’Time (seconds)’);ylabel(’Yaw Rate, r, (deg/sec)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED r RESPONSE’’);']J; 
eval (ans) 
text(.6,.85,’'* Measured Data Points ’,’sc’);pause; 
Peet rtcac*OUTZ (:,3),'og’ ); 
text(.6,.80,’o Estimated Response ’,’sc’);pause; 
if exist (’truth’)>0; 
Prowt, rede*~TRUZ(:,3),'—b’); 
fexti(.6,.750,'— “True Response" ‘,’sc’ ); 
end; pause; %meta A:\plots\OUTr 
% Bank Angle vs Time 
meee onrt;plot (t, rtdc*uydata(:,7),’*r’);hold on; 
xlabel(’Time (seconds)’);ylabel(’Bank Angle, phi, (deg)’); 
ans=[’/title(’’’,typac,’ ESTIMATED AND MEASURED phi RESPONSE’’);']; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’,’sc’);pause; 
meow, rLdc*OUTZ (:,4),'og’); 
text (.6,.80,’0o0 Estimated Response ’,’sc’) ;pause; 
me eexast (’ truth’ )>0; 
Prom e rede~rTRU2 (:,4) ,—b’ )7 
mexte(.6,.15,'—= “True Response" ’,'’sc’); 
end;pause; %*meta A:\plots\OUTphi 
%& Lateral G vs Time 
mempermort,;plot (t,uydata(:,8),’*r');hold on; 
xlabel(’Time (seconds)’);ylabel(’Lateral G, ay, (G)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED Lateral G 
RESPONSE’’);' Jj; 
eval (ans) 
text (.6,.85,’* Measured Data Points ’',’sc’) ;pause; 
mugen ,OUTZ(:,5),'’og’); 
text (.6,.80,’o0 Estimated Response ’,’sc’);pause; 
meexist (’ truth’ ) >0; 
fom (t, TRUZ(:,5),’—b’); 


mext(26,.75,’-— "True Response” ",'’sc’); 
end;pause; %meta A:\plots\OUTlatg 
moka Off ; 
%--------------------------------------- END MLEPLOT2 .M 
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E. ACTUAL LONGITUDINAL MMLE 


1. NPSMMLE3 .M 


% MACRO NAME >====== NPSMMLE3 .M ===S====< 
% Date: 31 Jan 91 
clear; 


ferase npsmmle3.1log; 
diary npsmmle3.log 


Gi Sp (0 srr re a ‘e; 
disp (’ NPS PARAMETER IDENTIFICATION MACRO FILE ') 

disp(’ FOR ACTUAL FLIGHT TESTS USING SIMPLFIED SHORT PERIOD ‘) 
disp (’ LONGITUDINAL STABILITY AND CONTROL DERIVATIVES’ ) 
disp(’ ’”) 

GAL SD (f mr 8 8 8 eo eee er 
npsinit3 %-------------------------------- RUN INITIALIZATION 
MACRO 


format compact, cle 
load npsinit3; 
global sref cbar gw iyy vtrue gqbar dt all ql thi anil del; 


§ me rn rrr rorerereee INPUT FLIGHT TEST DATA 
uydata=zeros (ndp, 6) ;% -------- ESTABLISH DATA MATRIX FOR MMLE.M 
$e rn rrr rrr nnn UYDATA(:,1) = DELTA E (INPUT) 
<a a ee ee oot UYDATA(:,3) = AOA 
eS UYDATA(:,4) = PITCH RATE (Q) 
%—-—---—--------------------------- UYDATA(:,5) = THETA 
$------------- -- --- -------------- UYDATA(:,6) = NORMAL ACC 
== — == COLUMN NUMBER ONE IN DATA FILE ELEVATOR INPUT 
coll=[’uydata(:,1)=simdata(:,1);’];eval(coll); 

$-—-—--—— COLUMN NUMBER TWO IN DATA FILE UNITY INPUT 
uydata(:,2)=ones (ndp,1); 

[SSS SS COLUMN NUMBER THREE IN DATA FILE ANGLE OF ATTACK 
col3=[’uydata(:,3)=simdata(:,2);']J;eval (col3) ; 

§----- COLUMN NUMBER FOUR IN DATA FILE PITCH RATE 
col4=[’uydata(:,4)=simdata(:,3);'’];eval(col4); 

5 COLUMN NUMBER FIVE IN DATA FILE THETA 
col5=[(’uydata(:,5)=simdata(:,4);'’J;eval(col5) ; 

$-—---— COLUMN NUMBER SIX IN DATA FILE NORMAL ACC 
col6=[’uydata(:,6)=simdata(:,5);'’]);eval(col6); 

%----SENSOR PLACEMENT CORRECTIONS FOR AOA AND NORMAL ACC 


uydata(:,3) =uydata(:,3)+(Xap*uydata(:,4) /vtrue) ; 
for 1=1idp, 

qsqr(1)=[uydata(1,4)]%2; 

Gz2-gqeqr’ ; 
end 
uydata(:,6) =uydata(:, 6) —(Zan*q2/32.17); 
qdot=[diff (uydata(:,4))*(1/dt) ;0]; 
uydata(:,6) =uydata(:, 6) —(Xan*qdot/32.17) ; 
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——. INITIAL CONDITIONS 
del=uydata(1,1);all=uydata (1,3) ;ql=uydata(1,4); 
thl=uydata (1,5) ;anl=uydata(l1, 6); 
i ADDITIONAL INPUTS TO MMLE FOLLOW 
peezeanam—'npsp2ss3’; %§ -------- MACRO NAME FOR P2SS FUNCTION 
p0O=pref; %- INITIAL PARAMETER ESTIMATES INPUT DURING NPSINIT3.M 
$-- CHECK IF THE WEIGHTING FUNCTION IS TO BE USED FOR INITIAL 
VALUES 
disp(’DO YOU WANT TO WEIGHT THE INITIAL ESTIMATES? INPUT 1=Y O=NO 
td 
) 
input (' '); 
if ans==1 
disp(’Input Weighting Row Vector length 1 x 5 ’) 
disp(’Use brackets- ex. [.1 11.11] & lower # higher weight’ ) 
rmsQ=input (’ ’); 


end 
pidgq=[1]);%--- IDENTIFY WHICH PARAMETERS ARE TO BE IDENTIFIED 
pidm=[1:5];%--- IN THE QUADRATIC, MARQUARDT, AND FINAL STAGES. 


pidf=[1:5];%- pidg, pidm, pidf MUST BE VALID EVEN IF NOT USED 
See=0 5 5 10 .02 .10 .001 1]; 

%- DEFAULT ITERATIONS AND CONVERGENCE CRITERIA, IF NOISE FREE 
OPT4=0 


ma —eve (4) *(.001) ; +------------- INNOVATIONS COVARIANCE MATRIX 
pert=le-4; $+-PERTURBATION USED FOR NUMERICAL GRADIENT CALCULATION 
linesearch=1;%--- USE LINESEARCH TO HELP PROC. NOISE CONVERGENCE 
els 

mmle $~-o-o-3 oo oro rrr rr CALL MAIN MMLE MACRO FROM TOOL BOX 
ooo PERFORM FINAL CALCULATIONS 


cla=pfin (1) ;cma=pfin (2) ;cmq=pfin (3) ; clde=pfin (4) ; cmde=pfin (5) ; 
cbicjs) ia) 

Gusp(’ ') 

deriv=({cla cma cmq clde cmde]; 

disp (’ MMLE STABILITY & CONTROL DERIVATIVES e) 
disp (’ CLA CMA CMQ CLDE CMDE’ ) 
disp (deriv) 

disp (’ INITIAL INPUT DERIVATIVES a 

disp (pref) 

pause 

mleplot3 

Gaary off 

$!print npsmmle3.1log; 

= a END NPSMMLE3 .M 
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2. NPSINIT3.M 


% MACRO FILE NAME >======= NPSINIT3.M =======< 
% Date: 31 Dec 91 
% INITIAL SETUP MACRO FOR RUNNING THE PARAMETER 
% ESTIMATION TOOLBOX IN MATLAB FOR SIMPLIFIED LONGITUDINAL 
% SHORT PERIOD STABILITY AND CONTROL DERIVATIVES 
% THIS MACRO ’GETS’ CONSTANTS AND DATA FILE 

$----- LOAD DATA FILE 

disp(’ DATA FILE MUST CONTAIN DATA IN COLUMN MATRIX '’) 
disp(’ MATRIX NAME SIMDATA: N DATA PTS X 5 COLUMNS ’ ) 

disp (’ ELEVATOR/ALPHA/THETA/PITCH RATE/NORMAL ACC.’ ) 


disp(’ *) 

aaisp DATA FILE NAME--MUST EXIST WITH A mat EXTENSION. os 
data=input (’ENTER DATA FILE NAME ( ==> WITH OUT <== .MAT 
EXTENSION)? ’,’8’); 


if exist (’ dt’ )==0,dt=input (’DELTA T BETWEEN DATA POINTS? ’);end 
ldc=[’load ’,data,’.mat;’']; 

eval (ldc) ;%--EXECUTES LOAD COMMAND 

ndp=length (simdata) ; 

t= (0 ndp=1 | ‘am: 

$------- INPUT REQUIRED CONSTANTS--- IF NOT IN DATA FILE 

if exist (’sref’)==0,sref=input (’REFERENCE AREA (S) IN SQUARE FEET? 
‘);end 

if exist (’ cbar’ ) ==0, cbar=input (’MEAN AERODYNAMIC CHORD IN FEET? 
’)y;end 

1f exist (’gw’ )==0,gw=input (’ AIRCRAFT GROSS WEIGHT IN POUNDS? ’);end 
if exist (/iyy’)==0,iyy=input (’MOMENT OF INERTIA (IYY) IN SLUG-FT*2? 
‘);end 

i1f exist (’ Xap’ ) ==0, Xap=input (’X-DIST FROM cg TO AOA PROBE (FT 
+FWD)’);end 

if exist ('Zan’)==0, Zan=input (’Z-DIST FROM cg TO NORMAL ACCEL (FT 
+DWN) ’) ;end 

if exist (’ Xan’) ==0, Xan=input (’X-DIST FROM cg TO NORMAL ACCEL (FT 
+FWD)’);end 

$sdc=[’save ’',data,’ simdata dt sref cbar gw iyy Xap Zan 

Xan’ ];eval(sdc); 

vtrue=input (’AIRSPEED IN FEET PER SECOND? '); 
altft=input (’ AIRCRAFT ALTITUDE IN FEET? oe i 

oat=input (‘OUTSIDE AIR TEMPERATURE? Ye 

‘---— se CALCULATE CONSTANTS DENSITY AND DYNAMIC PRESSURE 
rho=.0023769*exp ( (-1*32.174*altft) /(1716* (oat+460))); 
qbar=.5*rho*vtrue*vtrue; 
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—— INPUT INITIAL ESTIMATES FOR STABILITY 
Ko AND CONTROL DERIVATIVES 

pref (1)=input (‘CL ALPHA ESTIMATE (1/RAD)? '); 

pref (2)=input (‘CM ALPHA ESTIMATE (1/RAD)? '); 

pref (3)=input (’CM Q ESTIMATE (1/RAD)? '); 
pref (4)=input(’CL DE ESTIMATE (1/RAD)? '); 
pref (5)=input (‘CM DE ESTIMATE (1/RAD)? '); 
'erase npsinit3.mat; 

Save npsinit3 
a END NPSINIT3.M 


aie, 


3. NPSP2SS3.M 


function [a,phi, gam, c, 0d, aq 7ecae, Te nail bee SMe 


% MACRO FILE NAME >====>=>== NPSP2SS3. M =======< 

% Date: 31 Jan 92 

% cq cm cs sce cee gee GD ce ee ee ge ee ee oe ee ee ee ee ee ee ss oe a 
% MACRO TO EST. FUNCTION FOR TRANSFORMING 

% MODEL PARAMETERS INTO STATE SPACE EQUATIONS 

% a 
% P2SS FUNCTION FOR NPSMMLE3.M 

% a eae 
% p(l) = CL ALPHA | STABILITY AND CONTROL 

% p(2) = CM ALPHA | PARAMETERS 

% p(3) = CMQ | 

% p(4) = CL DE | 

% p(5) = CM DE | 

eee PERFORM INITIAL CALCULATIONS 


const1l=(-1*qbar*sref) / (cos (all) *gqw*vtrue/32.17) ; 
const2=qbar*sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 

const 4=32.17*cos (thl1) / (vtrue*cos (all) ); 
const5=(qbar*sref) /gw; 

const 6=-1* (const1*p(1) *all+ql+const1*p (4) *del+const4) ; 
const /=-1* (const2*p (2) *all+const3"p (3) *ql+eonst2*p (5) *dele 
const 8= (-1* (const5*p (1) *all+const5*p (4) *del))+anl; 


5 STABILITY DERIVATIVES 
a=[const1*p (1) il OF 

GOnstZ2*piZ) COnStS-p(3) > 

0 i OF lis 
NS SS SS SS SS SS SSS SS SS SSS CONTROL DERIVATIVES 
b=[const1*p (4) (const4+const6) ; 

COnStZ*~p (5) const7; 

0 (=) alii, 
a MEASUREMENT MATRIX 
c=[1 0 OL 

0 1 OF 

0 0 Ale 

COnses~pi (i) 0 Olin? 

a ee aa FEED THROUGH MATRIX 
d=[0 0; 

@) OF 

0 0; 

const5*p(4) COonses] ; 
ea ee oo STATE NOISE COVARIANCE 
q=eye (a) *le-4;%-- Q IS THE SAME SIZE AS a 
Sa eae et eel WITH Q*Q’ POS. DEFINITE 
ee ee ROWS IN Q IN WHICH PARAMETERS OCCUR, A VECTOR 
% SAME DIMENSION AS p 
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rowing=[0 0 0 0 0]; 


ee INITIAL STATE VECTOR 
xO0=[all ql thl]; 
0 DISCRETIZE 


foe “NEED TO EDIT dt (below) FOR THE DELTA T OF THE DATA ***x** 


dt=.1; 
[phi, gam]=c2d (a,b, dt) ; 
%------------------------------ END NPSP2SS3.M 


on: 


4. MLEPLOT3 .M 


% MACRO FILE NAME >=S=ss===== MLEPLOT3 .M =======< 
% Date: Sieevdan 92 
% =—--—---—- MACRO TO PLOT DATA FROM NPSMMLE3 


constl=(-1*qbar*sref) / (cos (all) *gw*vtrue/32.17) ; 
const2=qbar*sref*cbar/iyy; 

const3=const2*cbar/ (2*vtrue) ; 
const4=32.17*cos (thl) / (vtrue*cos (all) ); 

const5=(qbar*sref) /gw; 

const 6=-1* (const1*pfin (1) *all+ql+const1*pfin (4) *del+const4) ; 

const 7=-1* (const2*pfin (2) *all+const3*pfin (3) *ql+const2Z2*prinweeeeee 


const 8=(-1* (const5*pfin (1) *all+const5*pfin (4) *del) )+anl1; 


Se ae eee eee STABILITY DERIVATIVES 
a=[const1*pfin (1) 1 0; 

GOnstZ pLrm (2) const3*pfin (3) Gy. 

@) 1 Ol: 
ae ae ae oo CONTROL DERIVATIVES 
b=[const1*pfin (4) (const4+const6) ; 

COnStEZ prin, Conse. 

0 (-1*ql)]; 

: Sotueteetentnatententnatentnetententantentontententententnatnatnatententantententetentententententententontetan MEASUREMENT MATRIX 
c=[1 0 OF 

0 1 OF 

0 0 se 

CONStS pian) 0 Ole: 

Saha ee FEED THROUGH MATRIX 
d=[0 O- 

0 O- 

0 0; 

const5*pfin (4) const8]; 
rtdc=(180/pi) ; 

{=== == ———=--—— OUTZ = OUTPUT VECTOR---- OUT3 = STATE VECTOR 
[OUTZ, OUT3]=lsim(a,b,c, Gd, Uydata(*; 1-27 ee 
§=<<—-—------- PLOTS TO MONITOR AND STORED IN META FILE 


if exist ('typac’)==0,typac=input (’ INPUT THE AIRCRAFT TYPE ? 
Y,'s');end 

'erase a:\plots\outputg.met; 

'erase a:\plots\outputth.met; 

f'erase a:\plots\outputde.met; 

'erase a:\plots\outaoa.met; 

ferase a:\plots\outputq.met; 

$----ELEVATOR VS TIME 

hold off;ploeie, rede-uydatatiy yer )e- 

xlabel(’Time (seconds)’);ylabel(’Elevator Input (degrees) ’); 
ans=['title(’’’,typac,’ ELEVATOR INPUT VS TIME '’);']; 

eval (ans) 

pause 

meta A:\plots\OUTPUTde 
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$----- AOA (OBSERVED AND ESTIMATED) VS TIME 

Pemete(t, rtdac*uydata(:,3),’*r’);hold on; 

xlabel(’Time (seconds)’);ylabel(’AOA (degrees) ’); 
ans=['title(’’’,typac,’ ESTIMATED AND MEASURED AOA RESPONSE’’);']; 


eval (ans) 

text (.6,.85,’* Measured Data Points ’',’sc’);pause; 
Peete, rCAcC*OUTZ(:,1),'og'); 

text(.6,.80,’'o Estimated Response ’,’sc’);pause; 


pause 
meta A:\plots\outAOA 
pe === Q (OBSERVED AND ESTIMATED) VS TIME 


Mammon t, Plot (t, rtdc*uydata(:,4),’*r’);hold on; 

xlabel(’Time (seconds)’);ylabel(’Pitch Rate, Q, (deg/sec)’); 
ans=['’title(’’’,typac,’ ESTIMATED AND MEASURED q RESPONSE’ ’);']; 
eval (ans) 

text (.6,.85,’* Measured Data Points ’,’sc’);pause; 

eee, LCAC*OUTZ (:,2),'og'); 

text (.6,.80,’0 Estimated Response ’',’sc’);pause; 


pause 
meta A:\plots\OUTPUTQ 
$—-—---— THEATA (OBSERVED AND ESTIMATED) VS TIME 


Wemmemeor, - plot (t, rtdc*uydata(:,5),’'*r’') shold on; 

xlabel(’Time (seconds)’);ylabel(’Pitch Angle, Theta, (deg)’); 
ans=['’title(’’’,typac,’ ESTIMATED AND MEASURED THETA 
RESPONSE’ ’);'); 

eval (ans) 

text (.6,.85,’* Measured Data Points ’,’sc’);pause; 

pase, CCAC*OUTZ (:,3),'’0g'); 

text (.6,.80,’o Estimated Response ’,’sc’);pause; 


pause 
meta A:\plots\OUTPUTTH 
as ACCELERATION (OBSERVED AND ESTIMATED) VS TIME 


Memmeelorr;plot (t,uydata(:,6),’*r’);hold ony 

xlabel(’Time (seconds)’);ylabel(’Acceleration, G '’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED G RESPONSE’’);']; 
eval (ans) 

text (.6,.85,’* Measured Data Points ’,’sc’) ;pause; 

Peet, OUTZ(:,4),'o0G'); 

text (.6,.80,’o Estimated Response ',’sc’);pause; 

pause 

meta A:\plots\OUTPUTG 

hela® off; 

5S SS SS SSS SSS END MLEPLOT3 .M 
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FE. ACTUAL LATERAL-DIRECTIONAL MMLE 


1. NPSMMLE4 .M 


% MACRO NAME >=ss==== NPSMMLE4 .M ==s==s==< 
% Date: 5 Feb 92 
clear; 


l'erase npsmmle4.log; 


disp (ans) 

disp(’NPS STABILITY PARAMETER IDENTIFICATION MACRO FILE’ ) 
disp (’ FOR LATERAL—-DIRECTIONAL a 
disp (’ STABILITY AND CONTROL DERIVATIVES ry 
disp ie) 

disp (ans) 

npsinit4 Se ae RUN INITIALIZATION MACRO 


format compact,clc; 
load npsinit4; 
global sref bbar gw 1ixx 1ixz izz vtrue gqbar dt betal rolll yawl 
ranglel ayl dal drl; 
a NF III INPUT FLIGHT TEST DATA 
uydata=zeros (ndp, 8) ;% --- ESTABLISH DATA MATRIX FOR MMLE.M 
% UYDATA(:,1) = DELTA Aileron (INPUT) UYDATA(:,5) = ROLL RATE 
(Pp) 
% UYDATA(:,2) 
% UYDATA(:,3) 
(phi) 
% UYDATA(:, 4) 
(ay) 
%---- COLUMN NUMBER ONE IN DATA FILE AILERON INPUT 
coll=[’uydata(:,1)=simdata(:,1);’]J;eval(coll); 
%---- COLUMN NUMBER TWO IN DATA FILE RUDDER INPUT 
col2=[’uydata(:,2)=simdata(:,2);’];eval(col2); 
uydata(:,3)=ones(ndp,1);% UNITY INPUT 
5S COLUMN NUMBER FOUR IN DATA FILE BETA (beta) 
col4=[’uydata(:,4)=simdata(:,3);’]J;eval(col4); 
%---- COLUMN NUMBER FIVE IN DATA FILE ROLL RATE (p) 
col5=[’uydata(:,5)=simdata(:,4);’J;eval(col5); 
%---- COLUMN NUMBER SIX IN DATA FILE YAW RATE (r) 
col6é=[’uydata(:,6)=simdata(:,5);’J;eval(col6) ; 
%$---- COLUMN NUMBER SEVEN IN DATA FILE ROLL ANGLE (phi) 
col7=[’uydata(:,7)=simdata(:,6);'’];eval(col7); 
$---- COLUMN NUMBER EIGHT IN DATA FILE LATERAL G (ay) 
col8=[’ uydata(:,8)=simdata(:,7);’J;eval(col8) ; 
§—— = SENSOR PLACEMENT CORRECTIONS FOR beta AND ay 
uydata(:,4)=uydata(:,4)—-(Xbp/vtrue) *uydata (:, 6); 
rdot=([diff (uydata(:,6))*(1/dt) ;0]; 
pdot=[diff (uydata(:,5))*(1/dt) ;0]; 


DELTA Rudder (INPUT) UYDATA (:, 6) 
UNITY INPUT UYDATA (:,7) 


YAW RATE (r) 
ROLL ANGLE 


BETA (SIDE SLIP) UYDATA(:, 8) LATERAL G 
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uydata (:, 8) =uydata(:, 8) —(rdot*Xay/32.17) + (pdot*Zay/32.17) ; 
%----- INITIAL CONDITIONS FOR da, dr, BETA, p, r, phi, ay 
betal=uydata (1,4) ;rolll=uydata (1,5) ;yawl=uydata (1,6); 
ranglel=uydata (1,7) ;ayl=uydata (1,8); 

dal=uydata (1,1) ;drl=uydata (1,2); 


a — ADDITIONAL INPUTS TO MMLE FOLLOW 
p2snam=’npsp2ss4!’ ; $----- MACRO NAME FOR P2SS FUNCTION 
p0=pref;%---- INITIAL PARAMETER ESTIMATES 


%$--CHECK IF WEIGHTING FUNCTION DESIRED 
disp(’Do you want to WEIGHT the initial estimates? INPUT 1=YES 
O=NO’); 
ans=input (’ ‘’); 
if ans==1, 
disp (’ Input WEIGHTING ROW MATRIX: x. 12") 
mee Uo BRACKBITo— ex. [.leige1 21102 .1 11 =.11)'); 
disp(’NOTE *** The LOWER the # the HIGHER the WEIGHTING!’ ); 


mmsO=input (’ ’); 

end;clc; 
eeaa—| 1) ;+--------— IDENTIFY WHICH PARAMETERS ARE TO BE IDENTIFIED 
Pmeem—| 1:12) ; +--------- IN THE QUADRATIC, MARQUARDT, AND FINAL 
STAGES. 
waee=—| > 12) ;s—-————— pidq, pidm, pidf MUST BE VALID EVEN IF NOT USED 
opt=[0 5 5 10 .02 .10 .001 1]);%- DEFAULT ITERATIONS AND CONVERGENCE 
% CRITERIA, IF NOISE FREE OPT4=0 
aemeve (>) ~ (.01) ; 4—-—-———-—-———-—-—— INNOVATIONS COVARIANCE MATRIX 
pert=le-4;%-- PERTURBATION USED FOR NUMERICAL GRADIENT CALCULATION 
linesearch=1 ; %----- USE LINESEARCH TO HELP PROC. NOISE CONVERGENCE 
mm] e%$-ee rr errr rr CALL MAIN MMLE MACRO FROM TOOL BOX 
$n nnn nnn nnn PERFORM FINAL CALCULATIONS 


CYb=pfin (1) ;Clb=pfin (2) ; CNb=pfin (3) ; Clp=pfin (4) ; 
CNp=pfin (5) ;Clr=pfin (6) ;CNr=pfin(7) ;Clda=pfin (8) ; 
CNda=pfin (9) ;CY¥dr=pfin (10) ;Cldr=pfin (11) ;CNdr=pfin(12); 


deriv=[CYb Clb CNb Erp CNp Cia 

CNr Clda CNda CYar Cid CNdr]; 
disp (’ MMLE STABILITY & CONTROL DERIVATIVES -) 
disp(’CY b Cab CN b Cl p CN p Clee ) 
disp (deriv(l1,:)) 
disp(’CN r Clda CNda CYar Cres CNdr’ ) 
disp (deriv(2,:)) 
durep (’ INITIAL INPUT DERIVATIVES 4} 


Gmepipref(1:6)); 

Gumeepreft (7:12) ); 

pause;mleplot4 

diary off 

%$!print npsmmle4.log; 
a ee END NPSMMLE4 .M 
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2. NPSINIT4.M 


clear; 

% MACRO FILE NAME >======= NPSINIT4.M =======< 

% Date: 5 Feb 92 

*% INITIAL SETUP MACRO FOR RUNNING THE PARAMETER 

*% ESTIMATION TOOLBOX IN MATLAB FOR LATERAL—-DIRECTIONAL 

% STABILITY AND CONTROL DERIVATIVES 

% THIS MACRO 'GETS’ CONSTANTS AND DATA FILE 

disp(’ ') 

diiso(’ DATA FILE NAME~--MUST EXIST WITH A .mat EXTENSION. '’); 
disp (’ DATA FILE MUST CONTAIN DATA IN COLUMN MATRIX’ ) ; 

diva (’ MATRIX NAME SIMDATA: N (data pts) X 7 COLUMNS’); 

disp (’ AILERON/RUDDER/BETA/ROLL RATE/YAW RATE/ROLL ANGLE/LAT G’); 
data=input (’ENTER DATA FILE NAME ( ==> WITH OUT <== .MAT 


EXTENSION)? ',!'s’); 
ldc=[’load ’,data,’ .mat;’];eval (ldc) ; 
ndp=length (simdata) ; 


if exist (’dt’ )==0,dt=input (’ INPUT dt BETWEEN DATA POINTS. ‘y;end 
t=(0- Map ar, 

(eS INPUT REQUIRED CONSTANTS 

if exist(’sref’)==0,sref=input (’ REFERENCE AREA (S) IN SQUARE FEET? 
‘)y;end 

if exist (’bbar’ )==0,bbar=input (’WINGSPAN IN FEET? '’);end 

if exist (’ gw’ ) ==0, gw=input (’ AIRCRAFT GROSS WEIGHT IN POUNDS? ’);end 
if exist (’1ixx’ )==0,1xx=input (’MOMENT OF INERTIA (Ixx) IN SLUG-FT%2? 
’)y;end 

1f exist (’1ixz’) ==0,1xz=input (’MOMENT OF INERTIA (Ixz) IN SLUG-FT%2? 
')y;end 

if exist (’izz’)==0,izz=input (’MOMENT OF INERTIA (Izz) IN SLUG-FT*2? 
‘);end 

if exist (’ Xbp’ )==0, Xbp=input (’X-DIST FROM cg TO BETA PROBE (FT 
+FWD)’);end 


if exist (’ Zay’ )==0, Zay=input (’Z-DIST FROM cg TO LATERAL ACC. (FT 
+DWN) ’);end 


if exist (’ Xay’ )==0,Xay=input (’X-DIST FROM cg TO LATERAL ACC. (FT 
+FWD)’);end 

$sdc=[’save ',data,’ simdata dt sref bbar gw ixx i1xz i12z Xbp Zay 
Xay’]; 


$eval (sdc); 

vtrue=input (’AIRSPEED IN FEET PER SECOND? '’); 
altft=input (’AIRCRAFT ALTITUDE IN FEET? '’); 

oat=input (‘OUTSIDE AIR TEMPERATURE? a) 
rho=.002Z23763%expii(-i422. 174*alt ft) /(1716* (oat +460) )); 
qbar=.5*rho*vtrue*vtrue; 

5 EE ee INPUT INITIAL ESTIMATES FOR STABILITY 

5 AND CONTROL DERIVATIVES TO START MMLE PROGRAM 
pref (1)=-.6;% input (’CY beta FROM WIND TUNNEL (1/RAD)? '); 
pref (2)=-.15;% input (’Cl beta FROM WIND TUNNEL (1/RAD)? '); 
pref (3)=.20;% input (’CN_ beta FROM WIND TUNNEL (1/RAD)? ’); 
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eet (4)=—.35;% input (’Cl p FROM WIND TUNNEL (1/RAD) ? 


‘); 
peer (5)=-.05;% input (’CN p FROM WIND TUNNEL (1/RAD)? ’); 
pref (6)=.15;% input(’Cl r FROM WIND TUNNEL (1/RAD)? ’); 
pref (7)=-.2;% input(’CN r FROM WIND TUNNEL (1/RAD)? '); 
pref (8)=.05;% input(’Cl da § FROM WIND TUNNEL (1/RAD)? ’); 
pref (9)=-.001;% input(’CN da FROM WIND TUNNEL (1/RAD)? ’); 


pref (10)=.175;%input (‘CY dr FROM WIND TUNNEL (1/RAD)? ’); 
pref (11)=.02;%input(’Cl dr § FROM WIND TUNNEL (1/RAD)? ’); 
pref (12) =-.075;%input (‘CN dr FROM WIND TUNNEL (1/RAD)? '); 
lerase npsinit4.mat; 

Save npsinit4 

%---------------------- -- --- ---- - ----------- END NPSINIT4.M 
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3. NPSP2SS4.M 


function [a,phi,gam,c,d,q,x0,dt, rowing,b]=npsp2ss4 (p) 


*% MACRO FILE NAME >===S==== NPSP2SS4.M ====S===< 

% Date: 5 Feb 92 

% SS SS Sr Se SS SS eee 
% MACRO TO ESTABLISH FUNCTION FOR TRANSFORMING 

% MODEL PARAMETERS INTO STATE SPACE BQUATIONS 

% ee a SS Sn aS SS i 
% P2SS FUNCTION FOR NPSMMLE4.M 

% st ts St re 
% p(l) = Cie eea aa | p(7) = CN x 
%$ p(2) = Cl_ beta [STABILITY AND CONTROL | p(8) = Cl da 
6 p(3) = CN beta | Pr peo) = Seleea 
$ p(4) =Clp PARAMETERS | poy = Cyiaz 
% p(5) = CN | | p(l1l) = Cl dr 
Sp (6) fae ila | | p(lZ2) = Civaa 
$emern-- eee PERFORM INITIAL CALCULATIONS 


const1l=(qbar*sref) / (gw/32.17); 
const2=qbar*sref*bbar; const2a=const2*bbar/ (2*vtrue) ; 
const3=(qbar*sref) /gw; 


§----------------- INERTIAL MATRIX 
In=[1 0 O- 70; 
QO ixxs=1x2 0: 
0 —1 x7 Zz Oe 
0 0 GO Visi 
$n nee PLANT 
an=[const1*p(1)/vtrue 0 -1 (32.17/vtrue); 
const2*p(Z) const2a*p (4) const2a*p (6) OF 
Const 2*p (3) const2a*p (5) const2a*p (7) OF 
0 ah 0 OO]; 
% 
a=inv (In) *an; 
% 
bn=[0 const1*p(10) /vtrue OF 
const2*p (8) Genst2*pil1) OF 
const2*p (9) GONStZ-ptl2) 0; 
0 0 Oj; 
% 
b=inv (In) *bn; 
% 
c=[1l 0 0 0? 
0 al 0 OF 
0 0 a 0; 
0 0 0 ile 
COnNSES “=p (il) 0 0 Oj; 
% 
d=[0 0 betal; 
0 0 © 1 alee 
0 0 yawl; 
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0 0 ranglel; 


0 Conses*p (10) aylj; 
% 
a el STATE NOISE COVARIANCE 
q=eye (a) *le—-4; $---~------ Q IS THE SAME SIZE AS a 
% WITH Q*Q’ POS. DEFINITE! 
.——_-_-_------------ ROWS IN Q IN WHICH PARAMETERS OCCUR, A VECTOR 
% SAME DIMENSION AS p 
rowing=0*p; 
§ eee nnn nnn a INITIAL STATE VECTOR 
xO0=(betal rolll yawl ranglel]; 
arn nnn eee DISCRETIZE 


%x**X**NOTE***** CHANGE dt TO THE ACTUAL DATA VALUE **** 
dt=.05; 

[phi, gam] =c2d (a,b, dt) ; 

a ss a we ee END NPSP2SS4.M 
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4. MLEPLOT4 .M 


*% MACRO FILE NAME >======= MLEPLOT4.M =======< 

% Date: 5 Feb 92 

$e ee ee = MACRO TO PLOT DATA FROM NPSMMLE4 
$———------- GENERATE THE PREDICTED DATA 

§ = = = = = = === CALCULATE STABILITY AND CONTROL MATRIX 


mass=gw/32.17; rtdc=(180/pi) ; const1=(qbar*sref) /mass; 
const2=qbar*sref*bbar; const2a=const2*bbar/ (2*vtrue) ; 
const3=qbar*sref/gw; 
% 
In=[1 0 0 O; 

O aixx -1ixz 0; 

QO -ixz i2z2z O; 

0 0 0 ee 


an=[const1l*pfin(1)/vtrue 0 —1 32 .17/ Meer 
COnstZ ptm (2) const2a*pfin(4) constZ2a*pfin (6) OF 
CONnSt2Z prima) constZ2a*pfin(5) const2a*pfin (7) OF 
0 He 0 Ohilig. 
a=inv (In) *an; 
% 
bn=[0 const1l*pfin(10) /vtrue 0; 
const2*pfin (8) COnst2 prin |) O07 
const2*pfin (9) COnstZ*priumqi2) 0; 
0 0 0); 
b=inv (In) *bn; 
% 
c= (al 0 0 0; 
0 1 0 OF 
0 0 1 OF 
0 0 0 i 
Consts-prin (1) 40 0 Oi 
% 
d=[0 0 betal; 
0 0 row. 
0 0 yawl; 
0 0 ranglel; 
0 const3*pfin (10) ayl]; 
% OUT2 = OUTPUTS ---- OUT3 = STATE VECTOR 


[OUT2Z, OUT3]=lsim(4,b,¢,d, Uydata (5, 1:3)7t,-0)7 

% PLOTS FOR VIEWING ON MONITOR AND STORE TO META FILE 
!erase a:\plots\*.met 
subplot (Z1l); plot (t, rtde“uydaea(-, 1); 
xlabel(’Time (seconds)’);ylabel(’Aileron Input (degrees) ’); 
ans=[’title(’’’,typac,’ AILERON INPUT VS TIME ‘'’);’'J);eval (ans) 
subplot (212) plot (t, nede~uydacaien 2s, 
xlabel(’Time (seconds)’);ylabel(’Rudder Input (degrees)’); 
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ans=[‘’title(’’’,typac,’ RUDDER INPUT VS TIME ’'');'];eval (ans) 
pause; %*meta A:\plots\INPUTDAR 

% Beta vs Time 

mmmplot (111) plot (t, rtde*uydata(:,4),'*r’);hold on; 
xlabel(’Time (seconds)’);ylabel(’Beta (degrees) ’); 
ans=[‘title(’’’,typac,’ ESTIMATED AND MEASURED Beta 
RESPONSE’’);']; 

eval (ans) 

text (.6,.85,’'* Measured Data Points ’,’sc’) ;pause; 

pameeete , YCAc*OUTZ (:,1),'09g'); 

text (.6,.80,’o Estimated Response ’,’sc’);pause; 

pause; %meta A:\plots\outbeta 

*% Roll rate vs Time 

fweomeemont ;DlLot (t,rtdc*uydata(:,5),’*r’);hold on; 

xlabel(’Time (seconds)’);ylabel(’Roll Rate, p, (deg/sec)’); 
ans=[‘title(’’’,typac,’ ESTIMATED AND MEASURED p RESPONSE’’);']; 
eval (ans) 

text (.6,.85,'* Measured Data Points ’,’sc’);pause; 

ome, LCAC*OUTZ(:,2),'0g'); 

text (.6,.80,’o Estimated Response ’,’Sc’) ;pause; 
pause;%meta A:\plots\OUTP 

% Yaw Rate vs Time 

femmemort;, plot (t,rtdc*uydata(:,6),’*r’');hold on; 

xlabel(’Time (seconds)’);ylabel(’Yaw Rate, r, (deg/sec)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED r RESPONSB’’);']; 
eval (ans) 

text (.6,.85,’'* Measured Data Points ’,'’sc’) ;pause; 

Peete, TLAc*OUTZ (:,3),'’09g9'); 

text (.6,.80,’o0 Estimated Response ’,’sc’) ;pause; 

pause; %meta A:\plots\OUTr 

% Bank Angle vs Time 

moememort ;- plot (t,rtdc*uydata(:,7/),’*r’);hold on; 

xlabel(’Time (seconds)’);ylabel(’Bank Angle, phi, (deg)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED phi RESPONSE’’);']; 
eval (ans) 

text (.6,.85,’* Measured Data Points ’,’sc’);pause; 

pamemeee, TCAC*XOUTZ(:,4),'0g’ ); 

text (.6,.80,’o0 Estimated Response ’,’sc’);pause; 

pause; %meta A:\plots\OUTphi 
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% Lateral G vs Time 

hold off;plot (t, uydata(?:,8),  “r ) ;-neotamer. 

xlabel(’Time (seconds)’);ylabel(’Lateral G, ay, (G)’); 
ans=[’title(’’’,typac,’ ESTIMATED AND MEASURED Lateral G 
RESPONSE )) ; iz 

eval (ans) 

text (.6,.85,'* Measured Data Points ’,’sc’);pause; 

pilot (t, OUTZ(:,5);, caiman 

text (.6,.80,’0 Estimated Response ’,’sc’);pause; 

pause; %meta A:\plots\OUTlatg 


nold ofE - 
& — eee END MLEPLOT4 .M 
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APPENDIX B OUTPUT 


SIMULATED OUTPUT 
1. Longitudinal 
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Figure B-1 A-4D Elevator Input 
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Figure B-2 A-4D AOA Response 
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Figure B-3 A-4D Pitch Rate Response 
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A4D ESTIMATED AND MEASURED THETA RESPONSE 
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Figure B-4 A-4D Pitch Angle Response 


A4D ESTIMATED AND MEASURED G RESPONSE 
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Figure B-5 A-4D Normal Acceleration Response 
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A-4D SIMULATED LONGITUDINAL RESULTS 


jepltel p (pid) pref cramer 2fcramer 
1.0000 3.4568 ay, Sele 0090 0.0231 
2.0000 —Q. sea3 -0.8000 070007 0.0018 
3.0000 —3 4717 Sear CCUG Ove S11 0.0795 
4.0000 0.3441 0.8000 On0102 0.0260 
3) UNO a0) SONS: =i oU00 0.0020 OrCueZ 
MMLE STABILITY & CONTROL DERIVATIVES 
CLA CMA CMQ CLDE CMDE 
3.4568 —-0.36a3 —3 4 7 a7 0.3441 -0.4905 
INITIAL INPUT DERIVATIVES 
Be SO) 0G —0. 8000-50000 0.8000 -1.5000 
"TRUTH DERIVATIVES" USED TO GENERATE DATA 
3.4500 —-Omse00 —3._ 600 0.3600 -0.5000 
b. NAVION 
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Figure B-6 Navion Elevator Input 
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Figure B-7 Navion AOA Response 


NAVION ESTIMATED AND MEASURED q RESPONSE 


*» Measured Data Polnts 
o Estimoted Response 
— True Response 


Pitch Rate, Q, (deg/sec) 


Time (seconds) 





Figure B-8 Navion Pitch Rate Response 
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Figure B-9 Navion Pitch Angle Response 


NAVION ESTIMATED AND MEASURED G RESPONSE 
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Figure B-10 Navion Normal Acceleration Response 
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NAVION SIMULATED LONGITUDINAL RESULTS 


a p (pid) pref cramer 2fcramer insens 
1.0000 4.4567 5.0000 Oro 2 51: O07 £6 OHO 227 
2.0000 =O), 7097 =O S000 0.0064 0.0184 0.0046 
3.0000 =—8 .7238 > seeo0000 O21638 U2 7538 Ur ce35 
4.0000 oh, ShsOe S000 Cisse Balk Ge 0573 ORIOLES 
>. 00100 =). 8596 =1 2 00Go OO07 00307 0.0041 


MMLE STABILITY & CONTROL DERIVATIVES 


CLA CMA CMQ CLDE CMDE 
4.4567 -0.7097 -8 .7238 0.3300 -~0.8596 
INITIAL INPUT DERIVATIVES 
5.0000 -0.5000 -15.0000 0.5000 -~1.0000 
"TRUTH DERIVATIVES" USED TO GENERATE DATA 
4.4400 -0.6830 -9.9600 OR5550 —-0.9230 
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Figure B-1ll UAV Elevator Input 
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Figure B-12 UAV AOA Response 
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Figure B-13 UAV Pitch Rate Response 
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VAY ESTIMATED AND MEASURED THETA RESPONSE 
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Figure B-14 UAV Pitch Angle Response 
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Figure B-15 UAV Normal Acceleration Response 
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UAV SIMULATED LONGITUDINAL RESULTS 


pid p (pid) pref cramer 2fcramer insens 
.0000 5.1314 6.0000 0.0989 0.2913 0.0888 
.0000 0 . Soe =O" 6.00 C 0.0244 O:0WwaAlrs 0. Ole 
.0000 —-3:4725 2-16, 0000 LaZo28 a. Lae 0.6143 
.0000 0.0473 0.5000 0.0248 OOF: 0.0282 
.0000 -0. 2000 motel, 6).0119 0.0096 w= 0263 0.0085 


MMLE STABILITY & CONTROL DERIVATIVES 


CLA CMA CMQ CLDE CMDE 
Se 314 =O 29590 ~3 42> 0.0473 —0O 2285 


INITIAL INPUT DERIVATIVES 
6.0000 =. GN Oe OO O00 0. S000 mt |. 8)18/,8, 


"TRUTH DERIVATIVES" USED TO GENERATE DATA 
5 . OdaO'0 SOO O10 ee Ae) 10 0... O77aG +0 .3536 
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2. Lateral-Directional 
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Figure B-16 A-4D Aileron and Rudder Inputs 
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Figure B-17 A-4D Beta Response 
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Figure B-18 A-4D Yaw Rate Response 
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Figure B-19 A-4D Bank Angle Response 
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Figure B-20 A-4D Roll Rate Response 
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Figure B-21 A-4D Lateral Acceleration Response 
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Bb. NAVION 
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Figure B-22 Navion Aileron and Rudder Inputs 
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Figure B-23 Navion Beta Response 
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Figure B-24 Navion Yaw Rate Response 
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Figure B-25 Navion Bank Angle Response 
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Figure B-26 Navion Roll Rate Response 
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Figure B-27 Navion Lateral Acceleration Response 
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c. UAV 
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Figure B-28 UAV Aileron and Rudder Inputs 
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Figure B~-29 UAV Beta Response 
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Figure B-30 UAV Yaw Rate Response 
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Figure B-31 UAV Bank Angle Response 
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Figure B-32 UAV Roll Rate Response 
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Figure B-33 UAV Lateral Acceleration Response 
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Figure B-34 F-14A Elevator Input 
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Figure B-35 F-14A AOA Response 
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Figure B-36 F-14A Pitch Rate Response 
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Figure B-37 F-14A Pitch Angle Response 
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Figure B-38 F-14A Normal Acceleration Response 
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Figure B-39 T-37 Elevator Input 
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Figure B-40 T-37 AOA Response 
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Figure B-41 T-37 Pitch Rate Response 
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Figure B-42 T-37 Pitch Angle Response 
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Figure B-43 T-37 Normal Acceleration Response 


T-37 LONGITUDINAL FLIGHT TEST RESULTS 


pid p (pid) pref cramer 2fcramer insens 
1262.00 Te 12S 5.0000 0.0266 Oe. U'5 2 0.0206 
27.0000 =O 7 oae =O. 5000 OF ON 0.0282 0.0049 
3.000 0m -3Z modze -20.0000 0.5082 2.0071 0.2574 
4.0000 OE O08 023500 0. 038% GS 5 0. Og a 
5.0000 -1.8404 =i S00 C202 G 0.0749 0.0077 
MMLE STABILITY & CONTROL DERIVATIVES 
CLA CMA CMQ CLDE CMDE 
7s 125 —Oo267) *=SsZ2 a7 2 0.1068 -1.8404 
INITIAL INPUT DERIVATIVES 
50000 -0 .S0005=20°-0000 0. 3500 -17 31000 


Slag 8 0) 


INITIAL DISTRIBUTION LIST 
No. Copies 


Defense Technical Information Center Z 
Cameron Station 
Alexandria, Virginia 22304-6145 


Library, Code 52 2 
Naval Postgraduate School 
Monterey, California 93943-5100 


Chairman, Code AA 1 
Department of Aeronautics and Astronautics 

Naval Postgraduate School 

Monterey, California 93940-5000 


Professor Richard M. Howard, Code AA/Ho 4 
Department of Aeronautics and Astronautics 

Naval Postgraduate School 

Monterey, California 93943-5100 


Professor Louis V. Schmidt, Code AA/Sc i 
Department of Aeronautics and Astronautics 

Naval Postgraduate School 

Monterey, California 93944-5100 


Lcedr. Robert G. Graham, USN 1 
1169 Paramore Drive 
Virginia Beach, Virginia 23454 


carn. l.J. Barnes 1 
Unmanned Aerial Vehicles Joint Project 

Code PEO (CU) —-UD 

Washington, DC PUscl—-1014 


Mr. Richard J. Foch 1 
Ms. Peggy Toot 

Naval Research Lab, Code 5712 

4555 Overlook Ave. SW 

Washington, DC 20374 , 


Mr. Paul Heinold il 
Marine Corps UAV Project Office 

Intelligence Systems (C2IU) 

Marine Corps Research, Development, 

and Acquisition Command 

Quantico, Virginia 22134-5080 


URS ah 


10. 


Mr. Jim Murray 
Room, 1022 


NASA-Dryden Flight Research Facility 
Edwards AFB, California 93523 


2 




















WILE.  OINIVUA LID UNE 
NAVAL. P32 STGRADUATE SCHOO! 
MONTEREY CA 93943-5101 


GAYLORD $ 























‘ : 
‘ ‘ 
e . . 
a . ‘ + s 
al roe er . P 
= ee ‘ F : 
‘ te 
. ; ] ' ‘ 
’ berg . ; . 
{ a > ad o) ’ z e 
' : 4 . 
“ofe vane ate ; : 
re wy j ‘ @e fe 
: : 
4 "o@ 4 i, ‘ ! 
UL o : i, i 
Pr mae ’ 
' 
‘ t 
“ 









‘ 
ARLE LY 


cst ete a aie 3 2768 000180505 











“eleasdatag 6, 
ter r Beato ia 


t Ath 
a 


































’ 
« 
s ' 
e e i ' 
‘ > 3 fo ’ t ‘ 
‘ aa . * 
e e ¢ ‘ 
bs |e te) 4 ee 
: . ‘ew ‘4p 8 ta ts 
. 9 ‘ te . ‘ a oe : 
, sa te . ‘ 
’ 44 
. 
it 
THO test 
Le ee pa 
teh Pew ft 
: es pa. « Hay 
’ Foret of ofp g Bored 
Pl eP ot are i Dane Fi stele ene © ret a 
s4 far 









gtivas 


ote 





' . 
Par mew, 


alesse 
fete Faas 


t 
Le ee 
Hetty, 










1 Deore J 
he i *>s0t et gg 

pa Ty erase 
ee oe 9 

wit 

“otek er 
© Fat Ags ot 
ate 








































+f oye, 
eetie Behe Bd 

‘ ‘ig 
TNS fr). 
Fos 


a ae Be] 
Merete moan 
r¢ 









‘i. 
Fo@ateiog 

















ote 
Hath aga is eo, Danse, 
ot tadl ey ee ae ray 
Cd & vives vie ae 
ee ar 
Loy rs 
fo daties Otek . 
eal Mae civad 






































‘ fe 4 
+ 1 Dm Peteng 0 
Lee S Jud ria 4 
Malighs ot ay gis 4 
“@ . MES Ge Wr tan) cay iy 
ms ’ ta ‘ 
. ' g 
. a 4 
Le ol +29 é 
doherty Ca oe 
e SY we. aed Tidy 
Liter eten von? 
are 





B oevecpie es 
a! 

tiated 

iar to 




































bop 
Soe se ae} er) 
7 tet s' Paha tor 
ere Grah ei es cigar gon ae 
« t fea, CA ee | 
» + la hat ey ee iow, Le | ’ 
hale Sint ; : CL IURUA ae lace eter ne 
ry ae ae ee ig : ee . a aang ean art: 
Gabotla pie Suk ter ny, Sar Bh ove 
Peay fie oye S Le ACTS ST eT | 
Mot Lady LaPeer ree) 


oe 
4 bin sty ,! 
on 7 oe Le Oe 
wf toby, 









eae 
oe rbot & 












tO wate 
OP srpye 
° 2 


© tary 
tat 
fi 












omy 
4 fy, eta 
oh dalie IC Zul (er a 







Me 
fhe bat, 


Laat y 
FF se BS bs 
as 


he, Laer Ba ae 
had, 






















“edt ein vot 
Corte ow eee > 
FO pty 
| od 


a 8 44 
"aba Psy 
CAE 



















eye ne Cee 

fre OF Wont Bo bng 

tod tae gt 

. od EINE Aen Ye woe Cr 

wot sy, @ fel, 

a ea LG a ot Bt 

a , 
pte bey 

«Pas, 
















rate, 
gine 















rete £2 
Tree Te 1 
ves 8.8 


a 





a | 
teneg 


Paes 
tee 


eM A eee 


r= 


ee 
aie y’ , 












theta 
eT ERS So tess 











oh et te oy 


syte 
ebtegesaadad 5 '8 















| at 
Ving 
it ee TESS) 















Le Cee ee 
"+ e 
"SEV Seeteneers ay ne ae 
Sees rt ay Ws Nea rats 
TAM rmeetgn ie Hote 
° gies renee He 
ger ESS Be te 
Stuns 



















NEY ee, 
























































































































































. 
« * “ 
a4 F . ee . 
» “the es 40 ot 7 = 
eles ii a Sed): hag 
‘ COL org ‘ ‘oo. aa es heer aes 
TH ae +e yuu ‘ hon ave oh ae (sua 7 an =. i" 
teeta Ld Thay | “ ae eh . . . 6 4 w& Wa, ’ ‘ : . * . 
arith 'a ary ‘ a . + 7 P 7 ‘ +4 . + 
‘' ‘ ‘ . aris ‘ ay t+ the . La | ee ae | e ‘ e : wee Fi . : 
‘. = on 1 cr a Pry Srey ’ oo4 4 ot + pea bi ve ° 
: Cy edn on weee te ' 4 - a or ong : + ’ a 
ie Daa res aa 8 i toe , “ 
2 “pee Pa s 1 ‘ 
Maa Rae ' t ‘ Fs . 5 
: boda Lents eel ‘ . . 
wets . . is 
vrctys V'tratery . A i a . , 
é A bwty rier “at e . e ‘ 
ACEP Set eny Speedy ae atte erateet ce es a a 
4 oa abel 14 ee CSS rs A ‘. Fs 7 
bee Bytes Wyse oT . ‘ vs . rar 
Peary o Ad arn a . oe. 
Ma ae ees a pte. 
VF Sheps cys baal Malt, Guy a are , J ; 
TG ney, Fe WV gs eet ; Fhe age ; 
by 4°00 dese - . . 
to arate s pias ce ae “sey a os ' aie “SG . 
Las Sev, a, PAP stray ye Le weet: Waviw Nags, Seest artery |, ’ > @e 
wage ew “Ete cs yee “| Us » Te ose BER eB) sip pe apes ripe * . = 
pape els Sa Vig is Sere <a Ass & ytepage * Fi 
"a4" Aa iid Bah Sees Tr eres ee laiec ed oPla sys ge ae See e . r| - . 
WG a re tyte ee ielaeee ‘ FAA ee oe mele ae . . . 
“ : pees LT iran SFTEPatate v ge, oe tore * 
“wy oe 7 UMS sae ' ‘ ay. q 
=FM U4 the N 8 a Pane otter Th Rae 4 oe t 2 c 7 
Ona S a Lhd PSO gia a qe. . a 
eG soe sSetetage ee ry ' a r P 
tarvege fat enons op ' a4 . . 
Ferorne Fa" eho . : 
“VFO Wego TEsye SQlibetaasmer , ‘ ‘ fine é 
H ee ee oe i mat as 
etalerry : ter v : 
i: 3 ar Gye Fury, 
a eS REST “Hh gtey 
EVER eres eaey teat t 
‘5 e6-Ma on Ses Cote 
 Gs)rgy Atewey grin Perey = 
ems FAS Yr Fey roe ae, 
i bid) ae Lares Yas 
Vth Hn bey oaee, Se Cr 
revary peg MP8s Aba, 
Urey 










ere es 









ra 
Sa ee ary 
ahs Saeed Ce me 
urge iF mg 7) 4 















Gree RE 
Lh tah vo 
ALAS er 


oF 
“were 
wary 
















Ure ee 
ee 
reat RS es 
- ate ow, . 































































































































































































































































































































































’ ' 
. aes, | 
, < 
te .“ ’ 
we - + et. 
~S bey a, > > 
: a a ‘ ’ . 
. (we ' 3 
oa ‘ * Ce een . ‘ . 
SNEAK eb atetatartit gtr Sg fy tthe Pee, ; 
may roe thal eta y vers MG ON vas re oe ate ary woe tat 
Pat wa ake ® ike’ sft <4 ‘ +. SU ee ge “¢ " * 
. Idea Cy « az tea ba $< ‘ ‘ 
are wee tytyE Ueee tines “ea Muae ie ah Yee eer rie ae hie ee a "4 
. LWA UE ype : std au aL Wat ae bir | her © we oy + . : 
a ee Weve stares USragty Palys Peaiiteegacae tN Be 8 pee g » 4 , og 
S22 ity ere y § be shane gs mae, ' 
ry “4 PO we arate, «4 Ee 
Ty gy my aid Sie NCE or ee bherars tas, 
Oe VET Mim Bey RO ere ed y: Brass 7 : 
PRET oe ature BAe. ade SIRE ae t Th : 
Oa, org ays a as ones : : 
=f a 
cs 
hae 
yp ers wy ot Sto ee alle ook gS, hd 
i id rs ed Bd os Tal Le oy Lee a 4 
nAthadae Ts SPA Gent ere, ve a dity 
Yee re=, FET AFP od gone RIOT Uo tte Atay v ve aty sy Med ke . 
Mieka bel Tat 2 ten ET uoRresingm gig Wet rae at By 
PATE as omy Oe, “BSah o>, , Py Piare wea Ps 
“793, brary eg ers ate “Yat ors a Le ery . PSUINT Sood ony % 
Cee was ayy FS a ret as ore gripe gt Sota,» 
i A SYTW ae =) fey ity aie: 7 rae “4 
sin year ah iiteunasr sche 
fe UT : PS od oT Et SPA Nite oy eg ; 
oe ean tg 4 ue SyTb Ew galas eg” My aiyre Pane Arte “ a 
a tecns 4 ae. fea Utah 9p oo oa TL 
* ’ 
Bg DY oa? aria. % 
wks a's ota, aoa ae 
ty rte Uae Sythe agty,* re 2 * : 
A ee ry Mey J” Vie ts anaes? ae A , 
Wee fesasy noe 
stad ah) oa, TONE gg re 
aC Ys oy I - a 
BEA pepe coy "rule : 
Leureitg “alyen rey Bipte wae ; 
Big Pn ee, TaAS” vale ~ 6 
aly aps ohn eS . ‘ 
She ears day! "% 
eR La yea yw hiyey 
Plesuny te " Ber ary 7 
esate Na ey eras 
Ae Loe ett avay 
= stare ie J 
lee ebay; gee 1 E> see te 
3" tak 
cid 3 te Ven) 4 a. 
ee 
ir “het 
Wir 
Meteo) a RAY fe 
Bre 
Rarye ye Dee abe 
Pr bytary Ng? 





Sytimn a: 
Freeney GeAe 











+e 


